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An  alternative  approach  for  reducing  ground  level  concentrations 
of  sulfur  dioxide  emitted  from  electric  power  plants  is  presented.  The 
approach  involves  the  commitment  of  generating  units  using  an  extended 
version  of  the  method  of  production  costing.  It  is  shown  that  this  method, 
in  addition  to  providing  a  strategy  for  meeting  the  Federal  Ambient 
Air  Quality  Standards,  can  be  used  as  a  planning  model  for  fuel  purchases, 
plant  siting,  and  stack  height  determinat-ions.  It  is  easily  adapted 
to  any  size  electric  utility  system. 

The  historical  background  leading  to  the  promulgation  of  emission 
rate  and  air  quality  standards  is  first  discussed.  The  advantages  and 
disadvantages  of  several  of  the  more  prominent  existing  methods  are  also 
discussed.  Arguments  for  a  national  mix  of  alternative  strategies  are 
presented,  as  well  as  some  viewpoints  by  prominent  persons  on  the  issue 
of  continuous  or  noncontinuous  emission  controls. 

This  alternative  approach  involves  the  use  of  a  steady  state,  long 
term  (a  month  or  more)  mathematical  dispersion  model.  The  dispersion  model 
is  derived  in  a  novel  manner  from  the  steady  state  gradient  transport 
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equation  and  will  be  shown  to  be  the  Gaussian  dispersion  model  used 
extensively  for  predicting  pollutant  concentrations  from  point  sources 
(stacks).  The  statistical  approach  to  the  Gaussian  dispersion  model  is 
discussed  and  the  fundamental  relationship  connecting  the  two  approaches 
is  shown  to  be  only  a  special  case  of  a  more  general  connection  discovered 
and  presented  here. 

The  theory  of  production  costing  is  extended  to  include  monotonically 
increasing  or  decreasing  nonlinear  functions  of  the  electrical  output  of 
the  generator.  This  makes  it  possible  to  use  the  method  of  production 
costing  to  predict  long  term  estimates  of  the  sulfur  dioxide  emissions  as 
a  function  of  the  commitment  order  of  the  various  fossil  and  nonfossil 
electric  generating  units.  Larsen's  technique  for  relating  maximum  con- 
centrations to  average  concentrations  is  incorporated  into  the  approach 
in  order  to  extend  its  usefulness. 

A  digital  computer  technique  is  proposed  to  determine  the  optimal 
commitment  order  of  the  generating  units,  mix  of  fuels,  and  sulfur  contents. 
A  widely  used  model  for  predicting  average  concentrations  near  the  source 
is  shown  to  be  in  error  and  is  corrected  here.  A  simple  numerical  technique 
for  mathematically  constructing  the  "effective"  load  duration  curve 
is  presented  and  shown  to  involve  only  convolving  of  the  units--never 
deconvolving.  This  technique  avoids  many  of  the  computational  problems 
present  in  previous  models.  Finally,  the  application  of  the  computer 
program  to  a  realistic  size  electric  utility  system  is  presented  and 
discussed. 
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CHAPTER  1 
INTRODUCTION 

The  emissions  of  various  forms  of  sulfur  compounds  into  the  atmosphere 
have  been  a  concern  of  man  for  some  time  now,  although  it  has  only  been 
during  the  last  few  decades  that  it  has  received  so  much  public  attention. 
Hesketh  (1972)  reports 

As  early  as  1300,  a  royal  decree  was  issued  in  London 
prohibiting  the  use  of  low-grade  coal  for  heating  because  it 
created  excessive  smoke  and  soot.  The  only  known  case  of  capital 
punishment  because  of  an  air  pollution  violation  occurred 
in  the  13th  century  when  a  Londoner  violated  this  order. 
Sulfur  in  fuels  burns  to  sulfur  dioxide.  In  1600  sulfur 
dioxide  was  the  first  chemical  to  be  specifically  recognized 
as  an  air  pollutant.  However,  it  was  not  until  about  1940 
when  air  pollution,  as  such,  became  important,  (p.  1) 

Kellog  et  al .  (1972)  state 

1)  Man  is  now  contributing  about  one  half  as  much  as 
nature  to  the  total  atmospheric  burden  of  sulfur  compounds, 
but  by  A.D.  2000  he  will  be  contributing  about  as  much, 
and  in  the  Northern  Hemisphere  alone  he  will  be  more  than 
matching  nature. 

2)  In  industrialized  regions  he  is  overwhelming  natural 
processes,  and  the  removal  processes  are  slow  enough 
(several  days,  at  least)  so  that  the  increased  concentration 

is  marked  for  hundreds  to  thousands  of  kilometers  downwind,  (p.  595) 

Stoker  and  Seager  (1972)  have  determined  that  approximately  one  half  of 
the  total  atmospheric  burden  of  sulfur  compounds  results  from  the 
combustion  of  fossil  fuels  (primarily  coal  and  oil)  in  the  production 
of  electric  power.  Thus,  while  recognizing  that  there  are  other 
pollutants,  as  well  as  other  sources,  we  will  confine  ourselves  in  this 
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work  to  the  consideration  of  the  emissions  of  sulfur  dioxide  from 
electric  power  generating  plants.  (See  Chapter  3  for  a  discussion  of 
the  emissions  of  sulfur  dioxide  from  the  combustion  of  coal  or  oil 
containing  elemental  sulfur  or  sulfur  compounds.) 

1 . 1  Historical  Background  of  Clean  Air  Legislation 

In  1963  Congress  passed  the  first  Clean  Air  Act,  which  was  amended 
several  times  during  the  next  few  years.  The  1967  Clean  Air  Amendment 
required  that  the  Department  of  Health,  Education  and  Welfare  (HEW)  set 
air  quality  criteria.  The  first  such  criteria  were  published  in  February, 
1969  for  sulfur  dioxide  and  particulates.  The  criteria  contained  data, 
available  at  that  time,  on  the  adverse  health  effects  of  pollutants  and 
control  techniques  for  controlling  the  pollutants,  but  did  not  contain 
standards. 

On  December  2,  1970  the  Federal  Environmental  Protection  Agency  (EPA) 
was  established  by  an  act  of  Congress.  Among  its  many  functions  is  the 
supervision  of  air  pollution  control,  which  formerly  came  under  the  auspices 
of  the  Department  of  HEW  through  the  National  Air  Pollution  Control 
Administration  (NAPCA).  Some  weeks  later  the  91st  Congress  passed  the 
Clean  Air  Amendments  of  1970.  On  December  31,  1970  the  President  signed 
them  Into  law  as  Public  Law  91-604.  These  amendments  are  often  referred 
to  as  the  Clean  Air  Act  of  1970. 

One  of  the  many  provisions  in  the  act  calls  for  the  EPA  Administrator 
to  establish  Federal  Air  Quality  Standards,  and  on  April  30,  1971  these 
standards  were  issued  and  published  in  the  Federal  Register  (see  EPA 
(1971a)).  These  standards  included  both  primary  and  secondary  standards 
as  called  for  in  the  Act.  Section  109  (b)  (1)  specified  that  primary 
ambient  standards  should  be  selected  which,  allowing  for  adequate  margins 


of  safety,  would  protect  public  health.  Section  109  (b)  (2)  specified 
secondary  ambient  standards  to  protect  public  welfare  from  any  known 
or  anticipated  adverse  effects  from  pollutants.  Table  1  suirmarizes 
these  standards.  The  secondary  standards  for  the  annual  arithmetic  mean 
and  24  hour  average  were  adopted  but  subsequently  dropped. 

Air  quality  standards  are  not  enforceable,  in  that  the  air  cannot  be 
told  to  be  clean  with  any  assurance  that  the  air  will  somehow  do  so. 
Therefore,  emission  standards  are  required  which  are  enforceable.  The 
Clean  Air  Act  of  1970  provides  for  emission  limitations  to  be  set  by  the 
individual  states  through  State  Implementation  Plans  (EPA  (1971b)).  These 
plans  are  to  set  emission  limits  on  all  man-made  sources  such  that  the 
ambient  air  quality  standards  are  met  within  three  years  of  the  approval 
by  the  EPA  of  each  state  implementation  plan,  with  possible  extension  to 
mid-1977  through  state  initiative.  By  mid-1972  the  states  had  submitted 
their  implementation  plans. 

By  this  time  it  had  become  obvious  to  those  in  the  EPA  and  in  the 
utility  industry  that  there  simply  was  not  enough  low  sulfur  fuel  or 
other  control  technologies  (see  Chapter  2)  available  to  meet  and  maintain 
the  national  ambient  air  quality  standards  as  called  for  in  the  1970  Act. 
In  the  Fall  of  1972,  the  EPA  initiated  the  "clean  fuels  policy".  This 
policy  was  designed  to  encourage  states  to  delay  or  relax  their  sulfur 
dioxide  emissions  standards  in  areas  where  primary  ambient  standards  were 
not  being  exceeded,  and  to  give  priority  for  low  sulfur  fuel  to  plants 
which  were  threatening  the  primary  ambient  standards.  In  addition,  this 
policy  asked  that  the  states  promulgate  revised  implementation  plans  that 
did  not  contain  emission  limitations  more  stringent  than  required  to  meet 
ambient  standards  (see  EPA  (1974)).  According  to  the  EPA  (1975) 
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Table  1.  National  Primary  and  Secondary  Ambient 
Air  Quality  Standards. 


Primary  Standards  in 
Micrograms  per  Cubic  Meter 


Secondary  Standards  in 
Micrograms  per  Cubic  Meter 


Annual 
24  hour 
3  hour 


80 
365 


1300 
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"This  program  has  resulted  in  42  million  tons  of  coal  being  made  legally 
acceptable  through  revisions  in  the  sulfur  regulations  of  State 
Implementation  Plans."  (p.  15) 

In  the  Fall  of  1973,  the  much-publicized  Arab  oil  embargo  further 
complicated  the  fuel  supply  and  pressure  began  to  increase  from  the 
electric  utility  industry  to  amend  the  Clean  Air  Amendments  of  1970  once 
again.  On  March  13,  1974,  Rep.  Nelsen  (R-Minn)  introduced  a  bill 
(HR  13464)  to  amend  the  Clean  Air  Act.  This  bill  proposed  13  amendments 
be  adopted  by  Congress.  The  EPA  opposed  two  of  the  amendments  but  all  13 
were  transmitted  on  behalf  of  the  Administration  by  the  EPA  Administrator 
with  a  letter  explaining  the  differing  points  of  view  on  the  two  amendments, 
(We  will  mention  these  differing  points  of  view  in  a  discussion  of 
alternatives  in  Chapter  2.) 

One  significant  change  that  would  occur  if  these  amendments  were 
adopted,  would  be  the  extension  of  the  compliance  date  for  meeting  ambient 
air  quality  standards  to  1985.  Anticipating  that  Congress  would  hold 
hearings  during  1975  or  1976  on  these  proposed  amendments,  and  that  the 
electric  utility  industry  would  be  asked  to  respond  at  the  hearings, 
the  electric  utility  industry  formed  the  Clean  Air  Coordinating  Committee 
(CACC)  in  early  1975.  This  Committee,  which  consists  of  representatives 
from  the  principal  utility  industry  trade  associations  and  both  public  and 
investor-owned  electric  utilities,  has  been  responsible  for  coordinating 
the  major  efforts  of  the  utility  industry  to  respond  to  Congress  or  the 
regulatory  agencies  for  information  regarding  the  views  of  the  utility 


'Letter  dated  March  22,  1974  from  Russell  E.  Train,  EPA  Administrator. 
to  President  Richard  M.  Nixon,  in  forwarding  the  proposed  bill,  "The 
Clean  Air  Act  Amendments  of  1974." 
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Industry  on  the  proposed  amendments  (see  CACC  (1975)).  It  should  be 
noted  that  as  of  this  writing  (June,  1976)  the  Congressional  hearings 
have  not  been  held,  but  are  tentatively  scheduled  for  Spring  1977  (after 
the  Presidential  election). 

We  complete  this  section  with  a  quote  from  the  CACC  (1975)  that 
sunmarizes  their  feelings  (and  they  represent  the  electric  utility 
Industry)  on  the  state  of  existing  air  pollution  regulations: 

The  field  of  clean  air  regulation  is  in  a  state  of  flux. 
Congress  is  currently  reexamining  the  Clean  Air  Act  and 
the  Environmental  Protection  Agency's  implementation  of 
that  Act  to  date.  EPA  is  also  reassessing  its  own  policies 
under  the  Act.  The  Federal  Energy  Administration  and 
EPA  are  about  to  implement  the  Energy  Supply  and  Environmental 
Coordination  Act  (ESECA).  The  Federal  Power  Commission 
is  charged  with  certain  responsibilities  with  regard  to  the 
adequacy  and  costs  of  national  electricity  supply.  The 
Treasury  Department,  the  Department  of  Commerce  and  the 
Department  of  the  Interior  have  broad  responsibilities 
with  regard  to  control  of  mounting  inflation,  alleviation 
of  growing  unemployment  and  wiser  development  and  use  of 
national  resources.  The  Energy  Resources  Council  has  on- 
going responsibility  within  the  Executive  Branch  for 
coordinating  national  energy  policy  with  other  federal 
policies,  such  as  protection  of  the  human  environment,  the 
economic  well-being  of  our  citizens,  national  defense, 
the  balance  of  trade  and  foreign  policy,  (p.  1) 

1 .2  Statement  of  Purpose  and  Chapter  Outline 

It  is  the  purpose  of  this  dissertation  to  present  an  economical 
alternative  for  control  of  ground  level  concentrations  of  sulfur  dioxide 
from  electric  power  plants  using  the  method  of  production  costing.  In 
view  of  the  ambiguous  status  of  clear  air  regulations  discussed  above, 
it  is  not  surprising  that  there  are  numerous  alternative  approaches  to 
the  control  of  emissions  of  sulfur  dioxide  and/or  ground  level  concentrations 
of  sulfur  dioxide  from  electric  power  plants.  We  intend  to  present  a 
novel  approach  that  is  an  outgrowth  and  extension  of  the  work  of  Sullivan 
and  Hi  1  son  (1975). 
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Since  there  are  so  many  approaches  to  the  control  of  sulfur  dioxide 
from  electric  power  plants,  we  will  devote  Chapter  2  to  a  discussion  of 
the  relevant  features  of  the  inore  prominent  alternative  approaches  in 
existence  today.  In  Chapter  3  we  will  present  a  novel  solution  to  the 
steady  state  gradient  transport  equation  that  will  lead  to  the  classical 
Gaussian  dispersion  model.  In  Chapter  4  we  will  present  the  method  of 
production  costing  and  extend  its  theory  to  include  certain  generalized 
functions  of  the  electrical  output  of  an  electric  generator.  Also,  in 
Chapter  4,  we  will  present  a  statistical  method  due  to  Larsen  (1959)  that 
extends  the  usefulness  of  the  approach  taken  in  this  work.  Finally, 
in  Chapter  5,  we  will  apply  the  model  developed  in  Chapters  3  and  4 
to  a  realistic  size  electric  utility  system.  This  will  involve  a  computer 
program  which  is  also  discussed  in  Chapter  5.  A  statement  on  notation 
ends  this  chapter. 

The  notation  used  in  this  dissertation  is  discussed  as  it  is  first 
encountered.  This  is  appropriate  since  it  is  not  necessary  to  introduce 
any  nomenclature  that  is  not  standard  mathematical  notation.  Two  items 
worthy  of  mention  are   1)  vectors  are  always  underlined;  and  2)  a  bar 
over  a  letter  indicates  average  value,  e.g.,  u  is  read  u  bar  and  is  the 
mathematical  symbol  for  average  or  expected  value  of  u. 
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CHAPTER  2 
ALTERNATIVES  FOR  CONTROL  OF  SULFUR  DIOXIDE  EMISSIONS 

In  this  chapter  we  present  a  brief  discussion  of  some  of  the 
relevant  features  of  the  more  prominent  approaches  to  the  control  of 
sulfur  dioxide  emissions  from  electric  power  plants.  These  control 
approaches  can  be  conveniently  divided  into  two  broad  categories: 
continuous  emission  controls  and  noncontinuous  emission  controls  as  is 
often  done  in  the  literature  (see  PEDCo  (1975)  and  EPA  (1975)). 

Continuous  emission  controls  can  further  be  divided  into  precom- 
bustion  processes  and  postcombustion  processes.  Likewise,  noncontinuous 
emission  controls  can  be  divided  into  fuel-switching  or  load-switching 
or  a  combination  of  both.  In  Section  2.1  we  will  discuss  briefly  the 
continuous  emission  control  techniques.  Then,  in  Section  2.2  we  will 
discuss  the  noncontinuous  emission  control  approaches.  Finally,  in 
Section  2.3  we  will  present  arguments  for  an  interim  (up  to  1985)  mix  of 
both  continuous  and  noncontinuous  emission  controls  of  sulfur  dioxide  from 
electric  power  plants,  and  present  the  method  proposed  in  this  work. 

2 . 1  Continuous  Emission  Control 

Continuous  emission  control  techniques  limit  atmospheric  emissions  of 
pollutants  (sulfur  dioxide)  to  some  fixed  or  permanent  levels,  usually 
expressed  in  pounds  of  pollutant  emitted  per  million  BTU  of  heat  input. 
These  levels  are  to  be  consistent  with  achieving  State  Implementation  Plan 
emission  limitations  and  with  maintaining  ambient  air  quality  standards 
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under  the  worst  case  meteorological  conditions.  The  continuous  emission 
control  approaches  can  be  separated  into  precombustion  and  postcombustion 
methods.  We  will  consider  the  precombustion  processes  first.  We  will 
discuss  both  the  use  of  naturally  occurring  low  sulfur  fuels  and  the 
processing  of  high  sulfur  fuels  before  combustion. 

2.1.1  Low  Sulfur  Fuel 


Naturally  occurring  low  sulfur  fuels  used  for  electric  power  generation 
include  low  sulfur  coal,  low  sulfur  oil,  and  natural  gas.  Natural  gas  is 
in  short  supply  in  the  U.S.  and  has  been  given  the  lowest  supply  priority 
for  use  as  a  fuel  by  electric  power  plants  by  the  Federal  Power  Commission 
(PEDCo  (1975)).  The  1973  Arab  oil  embargo  followed  by  dramatic  increases 
in  crude  oil  prices,  has  produced  concern  about  dependence  on  foreign  oil. 
This  development  has  meant  that  except  for  some  power  plants  already  burning 
oil  along  the  east  and  west  coasts,  low  sulfur  oil  is  of  limited  usefulness 
as  a  clean  fuel.  Therefore,  we  only  consider  low  sulfur  coal  as  a  viable 
source  of  low  sulfur  fuel. 

There  are  large  reserves  of  low  sulfur  coal  in  the  Northern  Great 
Plains,  as  well  as  in  the  East.  The  reserves  in  the  Northern  Great  Plains 
are  located  far  from  the  major  eastern  markets  and  pose  transportation 
problems  and  added  costs.  The  low  sulfur  coal  reserves  in  the  East  are 
subject  to  strip  mining  legislation  since  much  of  it  is  located  on  steep 
slopes  (EPA  (1975)). 

The  use  of  low  sulfur  coal  has  some  environmental  problems.  For 
example,  most  electrostatic  precipitators  in  use  today  to  control  parti- 
culate emissions  depend  on  the  presence  of  sulfur  to  decrease  the 
resistivity  of  the  particles.  The  removal  of  this  sulfur  by  burning 
low  sulfur  coal  5  therefore,  decreases  the  effectiveness  of  the  precipitators. 
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sometimes  by   a  factor  of  ten  {Oglesby  and  Nichols  (1970)).  The  use  of 
coal  as  a  combustible  fuel  also  creates  fly  ash  disposal  problems  and  dust 
emissions  from  coal-handling  activities. 

There  are  derating  problems  associated  with  the  conversion  to  low 
sulfur  coal.  Oil-fired  plants  could  only  be  converted  if  they  were 
designed  originally  to  burn  coal  and  subsequently  converted  to  oil,  or 
if  they  were  designed  for  dual  oil/coal  firing  (PEDCo  (1975)).  In 
addition,  even  some  units  burning  high  sulfur  coal  would  experience 
deratings  in  the  range  of  10  to  30  percent  because  of  limited  pulverizing 
capability  (Perry  (1974)).  For  new  plants,  of  course,  these  deratings 
do  not  occur  since  the  boiler  would  be  designed  for  low  sulfur  coal. 

The  sulfur  content  limitation  for  new  plants  is  defined  by  the  New 
Source  Performance  Standards"^  (see  EPA  (1971c)).  For  aJM  other  facilities 
the  EPA  (1975)  states 

The  emission  limitation  is  determined  through  consideration 

of  the  following  factors:  (1)  dispersion  characteristics 

of  the  source  and  the  surrounding  area,  (2)  background 

concentration  of  the  pollutant  [sulfur  dioxide],  (3)  total 

emissions  rate  of  all  sources  in  the  area,  and  (4)  expected 

rate  of  growth  in  emissions.  The  applicable  sulfur  content 

limitation  is  based  on  allowable  emission  rates.  If  the 

determination  of  the  allowable  emission  rate  is  correct 

and  the  supply  of  coal  with  the  necessary  sulfur  content  is 

assumed,  this  alternative  is  better  than  99  percent 

effective  in  assuring  that  air  quality  goals  are  achieved,  (p.  21) 

In  addition  to  the  increased  cost  due  to  capital  conversion  of 
existing  plants,  and  the  increased  effective  cost  due  to  plant  deratings, 
there  is  an  increased  fuel  cost  for  the  low  sulfur  coal.  The  fuel  cost 
increase  is  generally  more  acceptable  to  electric  utilities,  since  these 
costs  can  be  passed  on  to  the  consumer  through  the  now  familiar  "fuel 


This  limitation  is  1.2  pounds  of  sulfur  dioxide  per  million  BTU  heat 
input. 
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adjustment  clause."  Whereas  capital  costs  can  only  be  recovered  through 
rate  base  changes  which  necessitate  rate  hearings. 

In  summary,  it  is  seen  that  low  sulfur  coal  is  a  viable  option  for 
meeting  emission  limitations  and  ambient  air  quality  standards.  Its 
principal  constraints  are  supply  and  technical  problems  associated  with 
using  them  in  certain  types  of  boilers  (PEDCo  (1975)). 

2.1.2  Fuel  Oil  Desulfurization 

Fuel  oil  desulfurization  is  only  one  of  several  approaches  we  will 
discuss  that  involve  preprocessing  of  the  fuel  before  combustion.  Hence, 
it  is  a  precombustion  process.  The  removal  of  sulfur  from  crude  oil  has 
been  demonstrated  at  many  refineries,  e.g.,  Phillips,  Humble,  Exxon, 
Chevron,  Gulf,  Standard,  and  others,  all  have  installations  that  are  in 
operation  today. 

The  desulfurization  process  is  capital  intensive,  requiring  on  the 
order  of  50  to  70  million  dollars  for  installation.  In  addition,  there 
are  significant  operating  costs,  and  the  process  requires  energy  that 
amounts  from  5  to  10  percent  of  the  energy  contained  in  the  oil  itself 
(Nelson  (1973)).  The  principal  drawback  to  the  use  of  low  sulfur  oil  as 
an  alternative  to  meeting  environmental  constraints  is  its  supply.  The 
shortage  of  capital  and  unpredictable  costs  of  crude  oil  have  severely 
limited  plans  for  new  installations  (EPA  (1974)).  For  this  reason  it  is 
probably  not  a  viable  option,  except  as  mentioned  before,  at  those  plants 
on  the  east  and  west  coast  that  are  already  using  imported  crude  oil. 

2.1.3  Coal  Desulfurization 

This  process  involves  physical  coal  cleaning  (sometimes  called  coal 
washing)  and  chemical  coal  cleaning  (sometimes  called  solvent  refining). 
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This  latter  cleaning  process  is  relatively  new,  and  in  fact,  the 
availability  is  forecast  to  be  in  the  late  1980's  (EPA  (1975)).  For  this 
reason  we  v/ill  only  discuss  the  physical  coal  cleaning  process. 

This  process  generally  involves  three  steps.  The  first  step  is  to 
separate  the  pulverized  coal  by  size.  The  second  step  involves  "cleaning" 
in  a  wet  medium  where  the  difference  in  specific  gravity  between  clean 
coal  and  impure  coal  is  used  to  separate  them.  Finally,  the  fine  cleaning 
involves  the  use  of  a  froth  flotation  in  which  the  separation  is  achieved 
by  floating  the  fine  coal  away  from  the  coal  containing  sulfur  and 
other  minerals  (Bodle  and  Vyas  (1974)). 

The  limiting  factor  in  this  coal  washing  process  is  that  only  the 
pyritic  sulfur  (mineral  compounds  containing  sulfur)  is  removed.  The 
organic  sulfur  (bound  in  the  coal  molecule)  is  not  affected.  Even  so, 
the  EPA  (1975)  estimates  that  about  14  percent  of  the  U.S.  coal  reserves 
can  be  cleaned  to  conform  to  the  New  Source  Performance  Standards.  In 
addition,  the  EPA  (1975)  estimates  another  55  percent  could  be  cleaned 
to  emit  less  than  4.0  pounds  of  sulfur  dioxide  per  million  BTU  heat  input. 
(The  Tennessee  Valley  Authority  (TVA)  has  a  plant  that  has  been  assigned 
an  emission  limit  of  5.2  pounds"^,  while  the  American  Electric  Power  System 
(AEP)  has  a  plant  with  an  emission  limit  of  6.0  pounds^''"  per  million  BTU 
heat  input.) 

Coal  used  in  the  preparation  of  coke  is  presently  cleaned  and  so 
the  technology  is  available.  Also  energy  requirements  and  costs  are 
modest  (PEDCo  (1975)).  A  significant  drawback  is  that  only  about  14  percent 


Personal  communication  with  E.  David  Daugherty,  iVA,  Chattanooga,  TN, 
Personal  commumcation  with  John  C.  Hoebel ,  AEP,  Mew  York,  M.Y. 
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of  the  U.S.  coal  reserves  could  be  cleaned  to  meet  the  Mew  Source 
Performance  Standards  {1.2  pounds),  and  so  coal  cleaning  is  at  best 
an  interim  approach  or  the  first  stage  of  some  other  approach.  Also, 
10  to  30  percent  of  the  coal  cleaned  by  weight  is  refuse  and  must  be 
disposed  of  in  landfills  which  involves  environmental  problems.  It  is, 
nevertheless,  a  viable  option  for  the  near  future. 

2.1.4  Coal  Gasification  and  Liquefaction 

Many  processes  for  converting  coal  to  a  (low  sulfur)  liquid  or  gas 
are  in  existence  or  under  development  today.  Yet,  this  approach  is  best 
considered  as  a  long  term  approach  and  not  of  much  value  to  the  electric 
utility  industry  in  the  immediate  future.  This  will  be  established  by 
citing  two  recent  reports.  The  report  by  PEDCo  (1975)  states 

High-sulfur  coals  are  also  converted  to  clean  fuel  by 
gasification  and  liquefaction.  The  systems  in  use  today, 
however,  are  generally  inefficient  and  expensive;  many 
have  considerable  potential  for  adverse  environmental 
impact.  Problems  of  low  efficiency  probably  can  be 
solved,  principally  in  operation  of  combined-cycle 
systems  now  under  intensive  development,  yery   high 
costs  of  installation  and  operation,  however,  coupled 
with  shortages  of  the  required  manpower  and  materials, 
could  preclude  competitive  commercialization  of  these 
processes  between  now  and  1990.  (pp.  xii-xiii) 

Further,  the  EPA  (1975)  reports 

While  coal  gasification  [and  liquefaction]  processes  are 
currently  in  commercial  use,  the  wide-scale  commercial 
application  of  coal  conversion  is  dependent  upon 
improvements  in  technology,  as  well  as  easing  of  con- 
straints on  the  construction  and  expansion  of  a  synthetic 
fuels  industry.  Technological  advances  that  increase 
gasifier  efficiencies,  ameliorate  environinental  impacts 
of  the  process,  and  decrease  costs  must  occur.  On  the 
construction  side,  it  is  estimated  that  a  full-scale 
plant  would  require  3  to  5  years  and  1.5  million 
man-hours  to  construct.  This  is  a  substantial  commitment 
of  manpower.  Additional  manpower  with  special  training 
is  required  to  operate  the  plant.  Therefore  it  is 
concluded  that  coal  conversion  processes  will  have  little 
impact  [on  electric  power-  systems]  between  now  and  1985.  (p.  23) 


»y"T  ■^.■¥>^watMrUiitif^iw  mirt 


14 


As  has  been  demonstrated  by  the  above  references,  the  conversion 
of  coal  to  a  liquid  or  gas  for  widespread  use  by  the  electric  utility 
industry  is  not  expected  before  about  1985  to  1990.  Thus  it  is  not  a 
viable  approach  for  the  near  future.  (As  is  explained  in  Section  2.3, 
the  alternative  approach  presented  in  this  work  is  a  viable  option  for 
the  interim  period  from  now  to  about  1985- ) 

2.1.5  Flue  Gas  Desulfurization 

This  is  the  first  of  two  techniques  we  will  discuss  that  involve  the 
removal  of  sulfur  after  the  combustion  process.  Hence,  flue  gas 
desulfurization  (FGD)  is  a  postcombustion  continuous  control  approach. 
This  is  the  approach  the  EPA  most  strongly  favors  (see  Section  2.3). 
At  the  present  time  there  are  about  a  dozen  different  FGD  processes  which 
have  been  demonstrated  in  full-scale  commercial  practice  (Princiotta  (1972)) 
Only  limestone  scrubbing  has  been  successful  for  extended  periods  of  time 
on  coal-fired  electric  power  plants  (EPRI  (1975)). 

The  limestone  scrubbing  process  is  a  nonregenerable  process,  in 
that  it  produces  a  waste  sludge  that  must  be  disposed  of  in  an  environ- 
mentally approved  manner--usuany  by  ponding  or  landfills.  For  disposal 
by  ponding,  a  100  megawatt  coal-fired  unit  would  require  0.5  acre  at  a 
depth  of  50  feet  each  year  the  FGD  unit  is  in  operation  (PEDCo  (1975)). 
.  PEDCo  (1975)  reports 

The  disposal  of  scrubber  sludges  entails  potential 
pollution  of  land  and  water.  Surface  waters  such  as 
rivers,  streams,  lakes,  and  ponds  can  be  contaminated 
by  leaching  and  percolation  of  sludge  liquor  into 
ground  water  through  soil  and  sludge  storage  areas. 
Large  areas  of  land  could  deteriorate  from  storage  of 
considerable  amounts  of  sludge  materials  that  typically 
contain  50  to  75  percent  water.  This  land  could  be 
made  useless  by  nonsettling  characteristics  of  the 
sludge,  (pp.  77-78) 


15 


In  Japan  scrubber  sludges  are  oxidized  to  form  fiber  gypsum  used  in  the 
production  of  vvall board;  and  if  the  economic  incentives  were  present, 
the  same  might  be  possible  in  the  U.S.  (EPRI  (1975)). 

Regenerable  processes  produce  a  marketable  by-product,  usually 
sulfur  or  sulfuric  acid.  The  most  promising  of  these  are  the  magnesium 
oxide  (Mag-Ox),  the  catalytic  oxidation  (Cat-Ox),  and  the  sodium  solution 
scrubbing  (Wellman-Lord)  processes. 

The  operations  of  the  FGD  system  will  require  from  three  to  six 
percent  of  a  power  plant's  total  (gross)  energy  input  (PEDCo  (1975)). 
Early  FGD  processes  were  unreliable  because  of  chemical  and  mechanical 
problems;  however,  more  recently  installed  systems  are  achieving  better 
than  90  percent  availability  (EPA  (1975)).  In  addition  to  sludge  disposal 
there  is  the  possibility  of  increased  local  concentrations  of  sulfur 
dioxide  near  the  plant  brought  about  by  decreased  buoyancy  of  the  plume. 
This  can  be  overcome  by  reheating  the  stack  gas  before  it  is  emitted. 

The  principal  objection  by  most  electric  utility  planners  has  been 
the  immense  capital  costs  associated  with  the  installation  of  a  system 
they  do  not  feel  is  a  proven  technology  (see  Section  2.3),  Nevertheless, 
the  EPA  (1974)  feels  "...[FGD]  represents  the  most  practical  medium  or 
long  term  solution  to  the  sulfur  oxides  problem  for  a  large  number  of 
coal-fired  power  plants."  (p.  9)  In  this  work  we  consider  FGD  systems 
as  a  viable  approach  to  controlling  sulfur  dioxide  emissions. 

2.1.6  Fluidized-bed  Combustion 

This  approach,  which  entails  a  combustion  modification,  is  one  of 
the  most  py^omising  approaches  for  using  fossil  fuels  in  the  production 
of  electric  energy.  It  offers  the  very   real  possibility  of  greater  thermal 
efficiencies  than  conventional  boilers.  The  thermal  efficiency  of  a 
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conventional,  boiler  with  FGD  (scrubbing)  is  about  37  percent,  and  no 
substantial  improvement  is  expected.  For  pressurized  fluidized-bed 
boilers,  thermal  efficiencies  of  38  percent  for  first-generation  units 
and  up  to  47  percent  for  second-generation  units  are  predicted  (PEDCo 
(1975)). 

In  addition,  fluidized-bed  boilers  should  be  able  to  burn  many  types 
of  fuels  and  combinations  of  fuels.  For  example,  all  ranks  of  coal, 
blends  of  oil,  gas,  rejects  from  coal-cleaning  plants,  and  even  municipal 
refuse  are  possible  fuels. 

The  use  of  a  sorbent,  such  as  limestone,  reacts  with  the  sulfur 
dioxide  formed  during  combustion,  and  if  the  sorbent  is  continuously 
withdrawn  (once-through  sorbent)  the  sulfur  is  effectively  removed.  Also 
proposed  are  systems  that  use  separate  vessels  to  extract  the  sulfur  from 
the  sorbent  and  then  reuse  the  sorbent  (sorbent-regeneration) .  This  latter 
process  eliminates  the  disposal  problem  of  the  once-through  sorbent.  The 
coal  fed  to  a  fluidized-bed  combustion  process  does  not  have  to  be  as 
fine  as  that  fed  to  a  pulverized-coal  boiler  so  particulates  (fly  ash) 
will  be  more  coarse  and,  therefore,  more  easily  collected. 

Cost  estimates  by  the  EPA  (1975)  indicate  that  pressurized  fluidized- 
bed  systems  should  have  capital  cost  savings  of  15  to  20  percent  and 
operating  cost  savings  of  10  to  15  percent  when  compared  to  conventional 
boilers  with  a  FGD  process.  For  later  generation  pressurized  systems  the 
operating  cost  savings  may  be  25  percent.  In  addition,  space  requirements 
are  comparable  to  existing  plants. 

The  EPA  (1974)  has  a  program  underway  to  demonstrate  a  pilot-plant 
pressurized  system  by  1982.  It  is  hoped  that  by  the  late  1980's  widespread 
use  of  these  systems  will  be  possible.  Thus,  as  promising  an  approach  as 
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it  is,  the  fluidized-bed  combustion  alternative  is  not  considered  in  this 
work  as  a  viable  option  for  the  near  term  (up  to  1985).  In  the  next 
section  vie  will  discuss  the  noncontinuous  emission  control  techniques. 

2.2  Noncontinuous  Emission  Control 


The  control  approaches  discussed  in  this  section  have  been  referred 
to  as:  intermittent  control  systems  (ICS),  supplementary  control  systems 
(SCS),  dynamic  emission  controls  (DEC),  and  various  other  names  that  imply 
noncontinuous  control.  A  noncontinuous  emission  control  system  has  been 
defined  by  the  EPA  (1973)  as  a  "...system  whereby  the  rate  of  emissions 
from  a  source  is  curtailed  when  meteorological  conditions  conducive  to 
high  ground-level  pollutant  concentrations  exist  or  are  anticipated." 
(p.  25697)  In  effect,  the  control  is  accomplished  by  reducing  emission 
levels  selectively  rather  than  continuously,  based  on  pollutant  dispersion 
characteristics  to  assure  that  ambient  air  quality  standards  are  not 
exceeded. 

The  Tennessee  Valley  Authority  (TVA)  calls  its  noncontinuous  emission 
control  system  the  sulfur  dioxide  emission  limitation  program  (SDEL).  SDEL 
is  the  most  comprehensive  SCS  in  operation  today  (PEDCo  (1975)).  It  is 
controversial  but  workable,  as  evidenced  by  the  successful  operation  at 
their  Paradise  power  plant  in  Drakesboro,  Ky.  since  1969  (Montgomery  et  al , 
(1975)).  The  SDEL  program  has  been  installed  at  nine  TVA  power  plants, 
which  have  been  divided  Into  class  1  and  class  2  type  programs.  The  class 
1  programs  are  less  complex  and  operate  in  an  open-loop  frame;  whereas, 
the  class  2  programs  are  operated  in  a  closed-loop  mode.  The  feedback 
is  via  a  telemetering  system  that  transmits  real-time  data  from  sulfur 
dioxide  monitors  located  strategically  throughout  the  geographical  area 
affected  by  the  sulfur  dioxide  emissions.  A  typical  class  2  program 
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involves  (see  Montgomery  et  al .  (1975)) 

1)  a  site-specific  pollutant  dispersion  model 

2)  12  sulfur  dioxide  ambient  monitors  (field) 

3)  sulfur  dioxide  emission  monitors  (stack) 

4)  at  least  one  instrumented  meteorological  tower 

5)  a  rawinsonde  system 

6)  several  minicomputers 

7)  a  telecommunication  system 

8)  an  environmental  data  station  to  be  operated  by  a 
staff  of  five 

Both  the  class  1  and  class  2  programs  are  supported  by  the  TVA 
meteorological  forecast  center  in  Muscle  Shoals,  Ala. 

The  noncontinuous  methods  for  reducing  emissions  involve  fuel-switching, 
load-switching,  a  combination  of  the  two,  and  tall  stacks.  Tall  stacks 
are  obviously  a  continuous  method  once  they  are  built.  However,  since 
they  do  not  reduce  emissions,  and  have  been  proposed  along  with  SCS,  they 
are  generally  thought  of  as  a  noncontinuous  emission  technique. 

2,2,1  Fuel -switching 

This  approach  involves  switching  to  a  low  sulfur  fuel  when  meteorological 
conditions  are  conducive  to  high  ground  level  concentrations  of  sulfur 
dioxide.  This  necessitates  the  storage  of  alternate  fuels  and  the  ability 
to  effect  the  switch.  All  of  the  necessary  components  for  fuel -switching 
are  currently  available  (EPA  (1975)).  Some  limitations  involve  the 
ability  to  switch  rapidly  enough  to  avoid  violating  the  ambient  standards, 
increased  particulate  emissions,  and  insignificant  decreases  in  the  total 
annual  emissions. 
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Some  positive  points  are  that  this  approach  becomes  effective  as 
soon  as  the  ability  to  switch  fuels  is  installed;  and  it  represents  an 
efficient  use  of  the  scarce  low  sulfur  fuels. 

2.2.2  Load-switching 

This  approach  is  similar  to  fuel-switching  except  that  when  conditions 
are  such  that  an  emission  reduction  is  required,  the  electrical  load  on 
the  particular  plant  is  switched  to  another  plant  in  the  interconnected 
electrical  transmission  system.  Again,  the  components  are  commercially 
available.  This  system  along  with  tall  stacks  is  the  primary  basis  for 
the  SDEL  program.  It  should  be  pointed  out  that  load-switching  has  long 
been  a  practice  of  the  electric  utility  industry  for  economically  loading 
of  generating  units  and,  also,  for  recovering  load  that  has  been  dropped 
due  to  forced  outages  of  generating  units. 

Similar  to  the  fuel-switching  approach,  drawbacks  to  this  approach 
include  questions  as  to  the  reliability  of  predicting  when  control  is 
needed.  Total  emissions  are  not  reduced  and,  in  fact  may  increase,  since 
the  alternative  plant  may  be  using  a  "dirtier"  fuel  or  may  be  more 
heavily  loaded. 

2.2.3  Tall  Stacks 

Any  process  that  regularly  generates  a  sizeable  plume  must  use  a 
stack  (chimney).  Two  criteria  involving  "good  engineering  practice" 
have  evolved: 

1)  The  stack  should  be  high  enough  to  allow  the  plume  to  escape 
from  the  wakes  created  in  the  lee  of  local  buildings  and 
topographic  features. 
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2)  The  exit  gas  velocity  should  be  high  enough  to  allow  the 
plume  to  escape  from  the  wake  of  the  stack  itself. 
The  first  criterion  is  usually  met  by  building  the  stack  2-1/2  times 
the  height  of  the  nearest  building  or  topographic  structure,  The  second 
criterion  is  usually  met  with  stack  gas  exit  velocities  of  1-1/2  times 
the  maximum  expected  normal  wind  velocities  at  the  top  of  the  stack. 
Thus  a  stack  exit  velocity  of  40  meters  per  sec  (m/sec)  would  be  required 
for  a  wind  speed  of  60  miles  per  hour  (mph)  (PEDCo  (1975)), 

The  EPA  (1974)  has  stated 

era's  "Tall  Stack  Policy"  encourages  the  use  of  stacks 
conforming  to  good  engineering  practice,  which  is  a  function 
of  the  individual  facility  configuration  and  local  terrain 
features.  In  general,  this  policy  results  in  stack  heights 
sufficiently  tall  to  minimize  ground  level  effects  caused  by 
aerodynamic  wakes,  eddies,  and  downwash,  and  those  caused 
by  high  winds  during  neutral  atmospheric  stability  conditions. 
In  some  cases  good  engineering  practice  requires  a 
relatively  tall  stack  to  overcome  adverse  terrain  features. 

However,  use  of  excessively  tall  stacks  in  an  attempt 
to  avoid  reducing  emissions  merely  results  in  dispersion 
of  sulfur  dioxide,  sulfates,  and  acid  aerosols  over  wide 
areas.  Their  use  as  a  substitute  for  permanent  [continuous] 
emission  controls,  in  addition  to  a  harmful  effect  on  health 
and  welfare,  would  be  inconsistent  with  the  ...[Clean  Air 
Act],  (pp.  14-15) 

The  use  of  tall  stacks  has  frequently  been  misunderstood.  For 
example,  in  the  July,  1975  EPRI  Progress  Report,  a  feature  article  by 
Yeager  (1975)  is  entitled  "Stacks  vs.  Scrubbers."  Stacks  are  a  necessary 
part  of  any  emission  control  system,  continuous  or  noncontinuous.  We 
will  view,  therefore,  the  use  of  a  tall  stack  as  a  method  to  be  used 
in  conjunction  with  (not  in  lieu  of)  some  other  method  of  control. 

It  should  be  mentioned  that  under  present  Federal  regulations 
(EPA  (1971c))  noncontinuous  emission  controls  are  applicable  only  to 
power  plants  for  which  construction  was  begun  prior  to  the  Issuance  of  the 
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New  Source  Performance  Standards  on  August  17,  1971.  With  a  doubling 
rate  of  ten  years,  about  one  half  of  the  power  plants  in  1981  will  still 
be  eligible  to  use  noncontinuous  emission  controls.  However,  this  approach 
at  best  must  be  considered  as  an  interim  (up  to  1985)  method. 

Two  of  the  drawbacks  common  to  all  noncontinuous  emission  control 
strategies  to  date  are 

1)  It  is  difficult  to  enforce  a  noncontinuous  emission  control 
system  in  a  multi -source  environment,  since  one  source  cannot 
control  the  emissions  of  another  source. 

2)  Insignificant  decreases  in  total  emissions  are  effected.  Thus, 
increased  levels  of  suspended  sulfates  and  "acid-rain"  problems 
are  possible. 

Two  of  the  advantages  common  to  noncontinuous  strategies  are 

1)  Noncontinuous  controls  are  relatively  inexpensive  when  compared 
to  continuous  strategies. 

2)  The  control  becomes  effective  as  soon  as  the  program  begins 
operation,  and  lead  times  for  initiation  of  the  control  are 
significantly  less  than  for  most  continuous  controls. 

We  have  discussed  some  general  aspects  of  noncontinuous  emission 
controls,  and  some  of  the  specific  requirements  of  the  SDEL  program  which 
is  being  used  by  TVA  (see  Montgomery  et  al .  (1975)).  Other  models  which 
have  been  proposed  are  Sullivan  (1972),  Sullivan  and  Hackett  (1973), 
Gent  and  Lamont  (1971),  and  Lamont  et  al .  (1975).  All  of  these  models 
suffer  from  the  same  kinds  of  limitations  mentioned  above.  In  the  next 
section  we  will  discuss  a  method  proposed  in  this  dissertation  that  is 
actually  a  hybrid  of  continuous  and  noncontinuous  emission  controls,  which 
Incorporates  some  of  the  better  features  of  both  control  approaches. 
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2 . 3  Arguments  for  a  National  Mix  of  Alternatives 

In  Section  2.2  we  discussed  the  noncontinuous  emission  control 
approaches.  All  of  the  models  mentioned  were  short  term  models.  That  is, 
their  primary  control  effort  is  directed  toward  the  3  hour  and  24  hour 
ambient  air  quality  standards.  Because  these  models  only  reduce  emissions 
during  adverse  meteorological  conditions,  their  principal  drawback  is 
that  they  do  not  significantly  reduce  anissions  over  the  long  term. 
A  second  limitation,  as  mentioned  before,  is  the  difficulty  of  enforcement 
in  a  multi-source  environment. 

The  model  we  are  proposing  in  this  work  avoids  the  principal  objection 
of  noncontinuous  emission  controls,  in  that  it  does  significantly  reduce 
long  term  emissions  of  sulfur  dioxide.  In  this  sense  it  is  a  continuous 
approach.  Since  it  also  involves  load-switching,  fuel-switching,  and 
tall  stacks,  it  has  elements  of  a  noncontinuous  approach.  Hence,  this 
approach  is  best  described  as  a  hybrid  approach.  We  prefer  to  call  it  an 
economical  alternative  approach  to  all  the  other  methods;  keeping  in  mind 
it  is  only  meant  to  be  an  interim  approach  valid  until  about  1985. 

Our  model  uses  a  dispersion  model  (Chapter  3)  to  predict  long  term 
air  quality,  based  on  long  term  emissions  determined,  a  priori,  through  a 
probabilistic  production  costing  mechanism  (Chapter  4).  That  is,  we 
simulate  a  power  system,  use  a  dispersion  model  to  predict  air  quality 
for  some  future  period  (usually  a  month),  and  through  the  method  of 
production  costing,  select  the  optimal  mix  of  load-switching  and  fuel- 
switching.  Tall  stacks  are  not  a  necessity,  but  are  desirable  as  discussed 
in  Section  2.2.  Before  applying  the  method  of  production  costing,  we 
must  extend  its  theory  to  include  certain  other  generalized  functions  of 
the  electrical  output  of  electric  generators.  This  is  done  in  Chapter  4. 
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The  use  of  the  extended  production  costing  method  allows  us  to 
select  the  order  of  committing  units  based  on  minimal  incremental 
additions  of  concentrations  of  sulfur  dioxide  at  ground  level.  In  effect, 
this  is  load-switching;  but  since  it  is  done  beforehand,  it  avoids 
most  of  the  drawbacks  of  load-switching  discussed  in  the  previous  section. 

Since  with  this  model  we  specify  the  fuels  and  their  maximum  sulfur 
contents  beforehand,  we  can  consider  fuel  supply  constraints,  boiler 
types,  and  fuel  costs  in  determining  the  optimal  mix  of  fuels  (fuel- 
switching).  As  is  discussed  in  Chapter  5,  the  model  lowers  the  sulfur 
contents  of  the  fuels,  where  necessary,  to  operate  the  power  system  in  an 
environmentally  acceptable  manner.  When  environmental  constraints  are 
met,  the  model  operates  in  an  economical  (fuel  costs  minimized)  mode. 

Since  the  model  determines  the  optimal  mix  of  fuels  (actually  sulfur 
contents  in  a  specified  type  of  fuel)  beforehand,  it  is  useful  as  a  planning 
technique  for  fuel  purchasing.  In  addition,  effects  on  the  economical 
or  environmental  costs  to  the  power  system  can  be  studied  relative  to 
fuel  changes,  plant  siting  (location),  stack  heights,  weather  changes, 
and  many  other  variables. 

Before  developing  the  model,  we  want  to  present  some  arguments  by 
prominent  persons  in  favor  of  a  mix  of  alternative  approaches  on  a 
national  scale  rather  than  selecting  just  one  or  two  approaches. 

In  a  letter  to  Secretary  Morton,  Chauncey  Starr,  President  of  ERPI 

+ 
states 


The  primary  standards  can  be  met  by  flue  gas 
desulfurization  devices  such  as  "scrubbers",  or 
alternatively,  by  an  intermittent  control  system 


t 
Letter  dated  March  26,  1975  from  Dr.  Chauncey  Starr,  President  of  the 
Electric  Power  Research  Institute  (EPRI)  to  the  Honorable  Rogers  C,B. 
Morton,  Secretary  of  the  Interior. 
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involving  a  mixture  of  fuel,  operating  schedules,  and 
tall  stacks,  all  adjusted  to  meet  the  day-to-day 
conditions  of  the  region.  Scrubbers  are  roughly  ten 
times  more  costly  than  intermittent  control  systems. 
As  noted  in  the  attached  memorandum,  scrubbers  do  have 
the  advantage  of  removing  sulfur  from  the  flue  gas  and 
thus  minimizing  uncertain,  long-term,  lov/-level  but 
cumulative  effects  of  sulfur  dispersal  on  the  regional 
environment.  For  this  reason,  EPA  has  preferred 
scrubbers.  In  either  case— scrubbers  or  intermittent 
controls--the  criteria  for  public  health  can  be  met. 

In  the  memorandum  referred  to  above  by  Dr.  Starr,  entitled, 
"Issues  and  Conclusions  on  the  Use  of  Intermittent  Control,"  it  is  stated 
in  the  conclusion 

Finally,  an  effective  SOx  control  strategy  must 
recognize  that  near-term  achievement  and  maintenance 
of  existing  SO2  ambient  air  quality  standards  depends 
on  understanding  the  following: 

a.)  Source  impact  on  ground  level  ambient  SO7  con- 
centration is  not  directly  related  to  mass  emission. 
This  encouraging  circumstance  provides  the  basis 
for  allocating  the  limited  quantity  of  low 
sulfur  fuel  and  control  technology  available  to 
the  specific  sources  where  it  is  most  needed  to 
achieve  and  maintain  the  standards. 

b.)  Intermittent  control  must  be  applied  to  the 

majority  of  large  utility  and  industrial  sources 
not  only  because  it  is  the  most  cost  effective 
approach  but  because  it  is  the  only  method 
which  can  be  made  available  in  sufficient 
quantity  in  the  near-term. 

c.)  Extensive  imposition  of  stack  gas  desulfurization 
for  controlling  sulfates  is  unjustified  not 
only  because  the  standard  has  not  been 
established,  but  because  the  current  state- 
of-the-art  is  technically  ineffective  for  this 
control  function,  (pp.  10-11) 

William  Lalor,  Senior  Vice  President  of  Southern  Services,  Inc.  in 
reference  to  Russell  Train,  the  EPA  Administrator,  says 


+ 

'Statement  of  William  G.  Lalor,  Jr.,  Senior  Vice  President,  Southern 

Services,  Inc.,  before  the  U.S.  House  of  Representatives,  Committee  on 

Interstate  and  Foreign  Commerce,  April  2,  1974. 
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Mr.  Train  seems  to  be  saying  that  intermittent  control 
systems  are  satisfactory  insofar  as  sulfur  oxides  are 
concerned  but  that  sulfur  oxides  are  not  the  problem. 
Therefore,  we  should  spend  a  great  deal  of  money  to  control 
sulfur  oxides  on  the  grounds  that  sulfates  at  some 
level  may  be  a  problem.  This  makes  no  sense  to  me. 
I  feel  that  Mr.  Train  has  set  up  a  sulfate  straw  man 
because  the  current  uniform  sulfur  oxide  emission 
standard  cannot  be  defended,  (p.  3) 

William  Donham  Crawford,  President  of  the  Edison  Electric  Institute, 

in  reference  to  a  meeting  with  Frank  Zarb,  Federal  Energy  Administration 

f 
Administrator,  states 

Mr.  Zarb  stated  that  the  Administration  did  not  wish 
to  repeat  last  spring's  situation  where  the  Administration 
had  recommended  several  amendments  to  the  Clean  Air  Act, 
but  EPA  refused  to  support  two  of  them  (intermittent 
control  systems  (ICS),  and  "no  significant  deterioration"). 
He  noted  that  the  proposals  had  languished  before  Congress 
ever  since  without  hearings  being  scheduled.  He 
thought  better  results  could  be  obtained  from  Congress 
if  the  Administration  presented  a  unified  front.  He 
said  tentative  agreement  had  been  reached  with 
Mr.  Train,  and  that  an  amendment  would  be  submitted 
in  January  along  the  following  lines:  Scrubbers  would 
be  designated  as  the  best  or  ultimate  technology 
and  permanent  emission  reduction  systems  should  be 
the  goal;  however,  it  would  be  recognized  that  more 
development  is  needed,  that  time  extensions  should 
be  allowed,  that  research  and  development  should  be 
continued,  that  ICS  should  be  authorized  in  the  interim, 
and  that  flexibility  should  be  permitted.  He  stated 
that  FEA  and  EPA  would  testify  in  support  of  the 
amendment,  and  that  the  industry  should  be  prepared  to 
present  its  views  fully  at  Congressional  hearings. 

The  meeting  above  was  held  in  late  December,  1974.  In  January,  1975 
the  Clean  Air  Coordinating  Committee  (CACC)  was  forFned  (see  Chapter  1). 

+4- 

George  C.  Freemen,  Jr.,  Special  Counsel  for  CACC  states  ' 


-f- 
Letter  dated  December  25,  1974  from  W.  Donham  Crawford,  President  of 
Edison  Electric  Institute,  to  Aubrey  Wagner,  Chairman  of  the  Board 
of  Directors  of  the  Tennessee  Valley  Authority. 

■  Letter  dated  April  30,  1975  to  the  Honorable  Frank  G.  Zarb,  Administrator, 
Federal  Energy  Administration,  from  George  C.  Freemen,  Jr.,  Special 
Counsel  for  the  Electric  Utility  Industry,  CACC. 
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Ue   believe  that  in  the  interim,  while  we  will  find 
the  long  term  answers  and  solutions,  the  government 
should  adopt  a  policy  that  requires  the  electric 
utility  industry  to  assure  attainment  of  present 
federal  air  quality  standards,  which  effectively 
protect  public  health  and  welfare,  but  at  the 
same  time  permits  each  plant  to  do  so  through 
whatever  mix  of  technology,  fuels  and  operating 
strategies  will  do  the  job.  This  would  minimize 
consumer  costs  and  contribute  to  this  nation's 
energy  self-sufficiency. 

Finally,  Yeager  (1975),  Program  Manager,  EPRI  states 

The  issue  of  tall  stacks  versus  scrubbers  is 
fundamentally  an  admini strati v.e  controversy  and  not 
a  technical  one.  This  is  the  main  reason  for  the 
delay  in  resolving  it.  In  reality,  a  technically 
feasible  national  strategy  for  achieving  the  esta- 
blished sulfur  dioxide  (SO2)  ambient  air  quality  standards 
over  the  next  decade  must  employ  both  stacks  and 
scrubbers.  Such  a  strategy  is  not  only  the  least 
costly  to  the  nation  but  the  only  one  consistent  with 
available  supplies  of  domestic  fuel  and  control 
technology,  (p.  2) 

In  the  same  article  Yeager  (1975)  states 

In  summary,  a  truly  feasible  national  SOn   control 
strategy  must  consider  the  roles  of  both  continuous 
and  intermittent  [noncontinuous]  controls  and  use 
established  techniques  to  apply  them  selectively.  Not 
only  will  such  a  strategy  effect  significant  cost 
savings  to  the  nation  but  it  is  the  only  one  that  can 
be  used  to  allocate  limited  clean  fuels  and  control 
technology  over  the  next  decade  without  sacrificing 
the  present  national  ambient  SO2  standards,  (p.  16) 

As  mentioned  above  (Crawford  meeting  with  Zarb)  the  EPA,  which  had 
opposed  noncontinuous  emission  controls  was  willing  to  change  its 
position.  We  will  close  this  chapter  with  statements  of  policy  from 
the  EPA  (1974)  and  EPA  (1975),  which  remain  in  effect  today  (June,  1976) 

EPA  (1974)  states 

...temporary  use  of  intermittent  controls  in 
carefully  selected  circumstances  will  facilitate 
more  rapid  attainment  of  the  current  primary  sulfur 
dioxide  standards  without  the  necessity  for  power 
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plant  shutdown,  will  allow  the  continued  usa  of  the 
nation's  high  sulfur  coal  reserves  while  control 
technology,  which  will  make  it  environmentally 
acceptable  is  being  installed,  and  will  also  allow 
time  to  increase  the  availability  of  low  sulfur  fuels. 
(p.  ii) 

Further  in  the  same  report,  the  EPA  (1974)  states 

Accordingly,  intermittent  control  is  currently 
considered  an  acceptable  control  measure  only  in 
cases  where  constant  emission  reduction  measures 
are  unavailable,  and  only  until  such  measures  become 
available.  Under  this  philosophy,  authorized 
intermittent  control  systems  are  referred  to  as 
"Supplementary  Control  Systems,"  meaning  that 
they  are  intended  to  supplement  available  constant 
emission  controls,  (p.  14) 

Finally  in  summarizing  the  report,  the  EPA  (1974)  states 

For  a  variety  of  reasons,  it  will  be  necessary 
to  apply  intermittent  controls  or  tall  stacks  tem- 
porarily on  some  plants  to  minimize  public  health 
impacts  from  sulfur  dioxide  until  adequate  emission 
control  measures  can  be  applied.  In  a  few  cases, 
interim  controls  may  be  required  until  the  early  1980's. 
The  use  of  intermittent  controls  will  be  minimized 
and  will  be  discontinued  as  soon  as  possible. 
Intermittent  controls  will  not  be  sanctioned  for 
long  term  use  where  constant  control  measures  are 
available,  (pp.  27-28) 

Then  in  1975  in  its  Report  to  Congress  on  Control  of  Sulfur  Oxides, 
required  by  the  Energy  Supply  and  Environmental  Coordination  Act  of 
1974. (ESECA),  the  EPA  stated  their  latest  position  (EPA  (1975)): 

EPA  believes,  hovjever,  that  a  number  of  isolated 
power  plants  can  use  supplementary  control  systems, 
which  require  additional  emission  reduction  on  an 
intermittent  basis,  as  an  interim  strategy  for  a 
number  of  years  without  increasing  the  risk  to 
public  health,  (p.  36) 

We  will  now  put  politics  aside.  In  the  next  chapter  we  begin  the 
theoretical  and  technical  development  of  the  alternative  model  we 
have  proposed. 


CHAPTER  3  , 

STEADY  STATE  MATHEMATICAL  DISPERSION  MODEL 

In  this  chapter  the  dispersion  model  used  throughout  this  dis- 
sertation is  developed.  While  it  would  be  desirable  to  have  a  dispersion 
model  capable  of  predicting  accurate  ground  level  concentrations  of 
sulfur  dioxide  emitted  from  electric  power  plants,  it  is  not  necessary 
that  the  estimates  be  exact.  This  is  because  decisions  about  the 
commitment  order  of  the  generating  units  (Chapter  4)  are  based  on 
comparing  alternative  strategies  (dispatching  orders)  and  selecting 
that  strategy  whose  incremental  increase  in  ground  level  concentrations 
throughout  a  geographical  area  is  the  smallest.  It  is  for  this  reason 
that  the  mathematical  dispersion  model  developed  in  this  chapter  need 
not  be  exact-- indeed,  the  physics  of  the  earth's  atmosphere  are  so 
complex  that  it  is  unlikely  that  an  exact  model  will  ever  be  found. 
If  the  dispersion  model  contains  a  bias  and  that  bias  is  consistent, 
then  the  bias  will  disappear  through  the  comparison  process.  Thus,  we 
are  concerned  with  the  effect  on  concentrations  of  one  strategy  relative 
to  another. 

A  gas  such  as  sulfur  dioxide  emitted  high  into  the  earth's  atmos- 
phere from  the  stack  (chimney)  of  an  electric  power  plant  will  undergo 


As  is  commonly  done  throughout  the  electric  power  industry,  we  will 
use  the  words  dispatching  and  committing  of  generating  units 
interchangeably. 
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considerable  dispersion  before  reaching  ground  level  (Turner  (1970)). 
This  dispersion  is  an  extremely  complex,  physical  mechanism  that  involves 
large  scale,  as  well  as  small  scale  turbulent  eddies  that  mix,  dilute, 
shear,  and  diffuse  the  gas  with  the  atmosphere.  Because  of  the 
complexity  of  turbulent  dispersion  in  the  atmosphere,  no  one  to  date 
has  proposed  a  mathematical  model  capable  of  describing  all  the  motions 
that  actually  take  place  (Tesche  et  al .  (1976)). 

There  are  two  general  approaches  to  developing  an  atmospheric 
dispersion  model.  The  two  approaches  are  the  gradient  transport  model 
credited  to  Pick  (1855)  and  the  statistical  model  formulated  by  Taylor 
(1921)  and  applied  by  Lin  (1960).  Both  approaches  have  had  wide  acceptance 
by  researchers  and  workers.  Both  have  some  limitations  due  mainly  to 
certain  simplifying  assumptions  in  the  interest  of  practicality.  Both,  as 
will  be  shown,  lead  to  a  Gaussian  dispersion  function  but  with  different 
second  moments. 

In  this  dissertation  the  gradient  transport  method  is  chosen  as  a 
starting  point  for  developing  the  dispersion  model.  We  shall  arrive 
at  the  classical  result  in  a  simple  but  novel  manner.  In  addition,  we 
shall  present  and  discuss  the  statistical  approach  and  show  it  is  identical 
to  our  derived  result  under  certain  conditions.  Finally,  we  shall  present 
the  plume  rise  formulas  used  with  the  dispersion  model  and  the  formula 
relating  gross  emissions  of  sulfur  dioxide  to  electrical  power  out  of 
the  generator.  Subsequently  in  Chapter  4  the  emission  formula  will  be 
modified  to  account  for  incremental  emissions  due  to  incremental  changes 
in  the  electrical  power  produced  by  a  generator. 
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3 • "!  Gradient  Transport  Model 

The  gradient  transport  model  is  chosen  because  it  can  be  derived 
from  basic  laws  of  physics.  What  happens  to  a  pollutant  after  it  is 
emitted  into  the  atmosphere  is,  indeed,  a  complex  question.  The  large 
scale  winds  move  it  downwind  from  its  source,  while  small  scale  winds 
(turbulent  eddies)  diffuse  and  mix  it  as  it  is  moved  downwind.  The 
atmospheric  temperature,  pressure,  and  humidity  influence  the  movement  as 
well  as  chemical  reactions  of  the  pollutant  with  other  agents  present. 
Sunlight,  clouds,  and  precipitation  cause  reactions  and  may  remove  some  of 
the  pollutant  from  the  atmosphere  or  convert  it  to  another  possibly  more 
dangerous  agent.  For  example,  the  conversion  of  sulfur  dioxide  to  acid- 
sulfates  is  considered  by  the  Environmental  Protection  Agency  (EPA)  to 
pose  a  serious  health  problem. 

A  precise  model  of  dispersion  would  have  to  include  all  of  the  above 
processes,  plus  other  weather  variables,  irregular  terrain,  variable  sources, 
and  types  of  deposition  (removal  from  the  atmosphere).  While  it  may  seem  a 
hopeless  task  at  this  point  to  develop  a  meaningful  model,  it  is  not.  This  is 
because  we  are  interested  in  concentrations  averaged  over  a  period  of  time. 
In  this  respect  we  are  fortunate  to  have  the  laws  of  statistics  in  our  favor. 
Although  the  individual  molecules  seem  to  move  in  a  completely  random 
fashion,  when  viewed  macroscopically  a  definite  plume  pattern  is  evident.  If 
one  recalls  viewing  a  visible  plume,  it  is  seen  to  flow  away  from  the  stack 
in  the  direction  of  the  transporting  winds  (large  scale),  while  it  spreads 
out  and  mixes  due  to  the  turbulence  of  the  eddy  wind  (small  scale).  This 
plume  pattern  that  emerges  is  due  to  the  laws  of  statistics  and  to  physical 
laws  such  as  the  conservation  of  mass,  momentum,  and  energy.  The  atmosphere 
must  obey  the  conservation  of  mass  for  both  the  pollutant  and  the  air. 
In  equation  form  (sometimes  called  the  equation  of  continuity) 
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this  is  written  as 


^$+  V  •  J  =  0  (3.1-1) 

dt 

where  x  is  the  concentration  of  the  pollutant  in  micrograms  per  cubic 
meter  (ug/m  )  and  J  is  the  flux  of  the  pollutant  in  micrograms  per 
square  meter  per  second  (ug/m  -sec). 

Just  as  ideal  gases  are  known  to  diffuse  from  regions  of  higher 
concentrations  to  regions  of  lower  concentrations,  the  gradient  transport 
theory  assumes  that  the  pollutant  flux  will  be  proportional  to  and  in  the 
direction  of  maximum  decrease  in  concentrations  per  unit  distance.  The 
proportionality  constant  is  called  the  diffusivity  constant  and  has  the 
units  square  meters  per  second  (m  /sec).  Stated  in  mathematical  terms, 
the  above  becomes,  for  one  dimension 

J  =  -K  1^  (3.1-2) 

X      X  9x 

where  K  is  the  proportionality  constant  (diffusivity  constant)  described 

A 

above  for  the  x-direction.  Equation  (3.1-2)  has  analogies  in  heat 
conduction,  neutron  flows,  and  molecular  motions.  It  is  credited  to 
Pick  (1855)  and  is  often  termed  Fickian  diffusion.  In  three  dimensions 
(cartesian)  equation  (3.1-2)  becomes 

J  =  -(^xl^ix^^|%^^zlfiz'  (3.1-3) 

where  (a^  ,  a  ,  a^  )  are  the  cartesian  system  unit  vectors  in  the  x,  y,  and 
X  — y   z 

z  directions , respectively.  For  an  isotropic  system  (K  =  K  =  K  =  K) 


"^Any  consistent  set  of  units  may  be  chosen.  Micrograms  per  cubic  meter 
has  been  chosen  by  the  U.S.  Government,  and  will  be  used  throughout  this 
dissertation. 
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equation  (3.1-3)  becomes  simply 

J  =  -KV  •  X  (3.1-4) 

Substitution  of  equation  (3.1-3)  into  (3.1-1)  yields 

It       3x  ^  X  8x  ^   9y  ^  y  dy   '        dz   ^  z  dz   '  U-l-b) 

Since  the  concentration  of  the  pollutant  may  change  in  time  as  well 
as  space  (due  to  wind  transport),  the  left  hand  side  of  equation  (3.1-5) 
may  be  expanded  as 


eft   9t   3x  dt   8y  dt   3z  dt 


The  concentration  at  a  particular  point  in  space  (x,  y,  z)  is  seen  to 
change  with  time  [dx/^t],   as  well  as  with  the  movement  of  the  pollutant 
with  velocities  dx/dt,  dy/dt,  and  dz/dt.  The  velocity  dx/dt  is 
the  velocity  with  which  the  pollutant  is  being  transported  in  the  x-direction, 
and  hence,  is  the  mean  wind  speed  u"  in  the  x-direction.  That  is 

-  dx 

U   =  -rr 
dt 

Similarly  v>fe  define 

^       dt 
the  mean  wind  speed  in  the  y-direction  and 

-  dz 
'^  =  dt 

the  mean  wind  speed  in  the  z~direction.  Defining  U  =  dx/dt  is  an 
approximation  because  dx/dt  is,  in  fact,  the  instantaneous  x-velocity 
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of  the  pollutant  which  would  include  the  small  scale  eddy  turbulence,  as 
well  as  the  large  scale  mean  wind  transport.  Again,  we  rely  on  the  law  of 
averaged  to  render  the  instantaneous  eddy  effects  insignificant  to  the 
mean  wind  effects,  since  in  the  large  they  tend  to  cancel  (Slade  (1968)). 
Their  effects  are  not  completely  ignored,  since  their  combined  effects 
are  included  on  the  right  hand  side  of  equation  (3.1-5). 

One  of  the  principal  limitations  of  &}]_  dispersion  models  is  the 
inability  to  describe  adequately  the  wind.  We  know  the  wind  changes 
speed  and  direction  in  time  and  space,  but  to  model  it  accurately  is  a 
problem  of  the  first  magnitude.  The  simplest  model  consists  of  describing 
the  wind  by  three  independent  components  averaged  over  time.  If  u 
represents  the  wind  speed,  it  may  be  written  as 

u_  =  U  a  +  7  a  +  w  a, 
X     ^'     — z 

where  u,  v,  and  w  are  the  components  of  mean  wind  speed  described  above. 
We  choose  to  use  this  rather  inaccurate  model  for  the  wind  because  we  are 
interested  in  concentrations  representative  of,  first,  steady  state,  and 
second,  long  term  (a  month  to  a  year).  As  a  result,  the  model  for  the  wind 
is  adequate,  and  in  fact,  is  the  one  used  by  most  workers  even  for  short 
term  concentrations. 

Below,  equation  (3.1-5)  is  rewritten  with  two  additional  terms,  R 
and  S  (Slade  (1968)) 

3t   "  3x  "^  ^  Sy  "^  ''  3z   3x  ^'^x  3x^   8y  ^\   9y ^  "  ^  ^^2  dz^^^^^ 

(3.1-6) 
The  presence  of  R  in  equation  (3.1-6)  is  to  include  the  rate  of  production 
(if  any)  of  the  pollutant  ciue  to  chemical  reactions,  S  is  the  source  of 
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emissions  (positive)  and  may  include  deposition  (negative)  at  the  position 
(x,  y,  z)  in  the  atmosphere. 

Equation  (3.1-6)  is  very   general  and  many  simplifications  are  needed 
to  solve  it.  It  is  the  basis  for  all  gradient  transport  studies.  It 
was  first  investigated  by  Schmidt  (1925)  and  Richardson  (1925).  Because 
of  the  K  terms  in  the  equation,  gradient  transport  theory  is  often  termed 
K-theory.  In  the  next  section  we  will  attempt  to  solve  equation  (3.1-6) 
with  appropriate  boundary  conditions. 

3.2  Steady  State  Gradient  Transport  Model 

It  is  instructive  to  consider  each  term  of  equation  (3.1-6).  Since 
in  this  work  we  are  concerned  with  the  steady  state  solution  to  equation 
(3.1-6),  we  immediately  set  the  first  term  on  the  left,  9x/9t  =  0.  The 
next  three  terms  on  the  left  are  the  transport  terms  due  to  large  scale 
mean  winds  (sometimes  called  advection  terms).  If  one  recalls  viewing 
a  visible  plume,  the  pattern  that  evolves  is  one  of  a  plume  streaming 
off  in  one  predominant  direction.  We  will  take  that  downwind  direction 
to  be  the  x-direction  (a^).  It  is  assumed  then,  that  the  horizontal  (a^  ) 
mean  wind  speed,  v,  and  the  vertical  (a_  )  mean  wind  speed,  w,  are  small 
compared  to  the  downwind  mean  wind  speed,  u.  Thus,  the  mathematical  model 
for  the  steady  state  wind  is  further  simplified  to 

u=ua^  (3.2-1) 

and  the  direction  is  taken  to  coincide  with  the  observable  downwind 
direction. 

Inspection  of  equation  (3,1-6)  indicates  there  are  two  dispersion 
effects  taking  place  in  the  x-direction.  The  first,  on  the  left-hand 
side,  is  the  advection  or  transport  term  due  to  the  mean  wind  speed,  u, 
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and  the  second,  on  the  right-hand  side,  is  the  diffusion  term  (K^-terni) 
due  to  the  turbulent  eddies.  In  the  downwind  direction,  the  advection 
term  far  outweighs  the  diffusion  term  under  normal  circumstances.  We 
assume  this  to  be  true  for  our  steady  state  model.  Mathematically,  this 
is  equivalent  to 

"If^^lri^li'  .  (3-2-2) 

In  the  horizontal  and  vertical  directions,  just  the  opposite  is  true. 
The  diffusion  terms  outweigh  the  advection  terms  significantly.  Thus 

7^«—  (K  ■^)  f3  2-3^ 

3y   3y  ^  y  dy'  ^-^'^  ^' 


and 


^!^«^(^l^)  (3.2-4) 


92   9z  '  z  3z 

Equations  (3.2-3)  and  (3.2-4)  are  consistent  with  the  assumption  leading 
to  equation  (3.2-1). 

The  R  term  in  equation  (3.1-6)  represents  the  chemical  reactions 
that  may  take  place  to  generate  or  remove  the  pollutant.  In  more  com- 
plicated studies  where  many  pollutants  are  being  modeled,  this  term  serves 
to  couple  many  equations  of  the  form  of  equation  (3.1-6),  thereby  greatly 
complicating  the  problem.  In  this  work  we  are  modeling  only  sulfur 
dioxide.  We  will  take  R  to  be  zero.  As  regards  the  removal  of  sulfur 
dioxide,  this  amounts  to  slight  overestimates,  which  means  being  conservative 
in  our  decision-making  process.  The  generation  of  sulfur  dioxide,  whether 
chemical  or  otherwise,  will  be  included  in  the  S  (source)  term. 
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The  final  term  in  equation  (3.1-6)  is  the  source  term,  S.  This 
term  represents  the  production  and  deposition  of  sulfur  dioxide  at  any 
point  (x,  y,  z),  in  the  atmosphere.  As  explained  above,  the  deposition 
is  ignored,  which  produces  conservative  estimates.  Electric  power  plant 
stacks  are  the  only  sources  of  sulfur  dioxide  considered  here.  This  is 
because  we  want  to  evaluate  and  compare  the  effect  each  plant  has  on 
adding  sulfur  dioxide  at  ground  level  throughout  some  appropriate  geogra- 
phical area.  Background  contributions,  such  as  area  sources  of  sulfur 
dioxide  can  be  included  afterwards  as  a  percentage  of  the  total  concen- 
trations. 

Given  a  system  of  n  stacks,  there  would  be  sources  at  the  following 
points  in  space 

(x-j,  y^,  z^),  (X2,  y2,  Z2),...,(x^,  y^,  z^) 

where  (x.,  y.,   2.)   is  the  location  of  the  i-th  stack  (i  =  l,2,...,n). 
Everywhere  else  S  would  be  zero.  The  simplest  method  of  handling  these 
multiple  sources  is  to  model  each  separately  and  then  superimpose  the 
concentrations  due  to  each  source.  That  superposition  is  valid  follows 
from  the  fact  that  pollutants  regardless  of  their  source  (location)  must 
follow  the  same  basic  dispersion  laws. 

The  use  of  the  superposition  principle  means  equation  (3.1-6)  need 
only  be  solved  once  for  a  given  location  of  the  cartesian  coordinate 
system.  The  solution  for  another  source  at  a  different  location  simply 
requires  a  coordinate  system  translation  and  a  new  source  strength. 

For  a  single  source  located  at  the  origin  (0,0,0),  and  all  of  the 
above  simplifications,  equation  (3.1-6)  reduces  to 
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u^ 


3x   ar^l^'^'^^if^^^^Q'Q'^'^)  (3.2-5) 

where  S(Q, 0,0,0)  is  the  source  of  strength,  Q.  (yg/sec)  located  at  the 
origin.  Instead  of  including  the  source  in  the  dispersion  equation,  as 
it  is  in  equation  (3,2-5),  it  may  be  included  in  a  boundary  condition. 
The  appropriate  boundary  condition  requires  that  the  concentration  at  any 
point  in  space  moving  (at  speed  u")  across  an  infinite"^  zy-plane  must  equal 
the  source  strength,  Q.  Mathematically  this  is  equivalent  to 


»  00  .  00 


ux  dydz  =  Q  •  (3.2-6) 

-ooJ  -oo 

and  is  another  statement  of  the  conservation  of  mass.  There  are  two 
other  boundary  conditions  that  are  necessary.  The  concentrations  must 
approach  zero  at  very   large  distances  from  the  source,  and,  the  concen- 
tration must  be  infinite  at  the  source. 

The  remaining  assumption,  before  actually  solving  equation  (3.2-5), 
involves  the  two  eddy  diffusivity  constants,  K  and  K  .  We  will  assume 
K^  to  be  independent  of  y  and  K  to  be  independent  of  z.  This  is  not 
tantamount  to  assuming  an  isotropic  atmosphere,  as  is  frequently  done. 
In  fact,  as  will  be  shown  later,  K  and  K  are  both  functions  of  the 
downwind  spatial  coordinate,  x.  Equation  (3.2-5)  and  the  three  boundary 
conditions  are  rewritten  below  with  the  above  assumptions  incorporated 

"|J-V4*^4  (3.2-7) 

ay  dZ 


t 


Physically,  the  yz-plane  would  have  to  terminate  at  the  earth's  surface 
Mathematically,  it  is  convenient  to  ignore  the  earth's  surface  in  equation 
(3.2-5),  and  to  include  its  effect  through  boundarv  reflections.  This  is 
discussed  in  Section  3.4. 
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u  X  dydz  =  Q  (3.2-7a) 


lim  X  =  0  (3.2-7b) 

r  ^  CO 


lim  X  =  -  (3.2-7c) 

r  ^  0 

2    2    2  1/2 
where  r  =  (x  +  y  +  z  )  '  .  Equation  (3.2-7)  and  its  associated  boundary 

conditions  will  now  be  solved. 

There  are  many  approaches  to  the  solution  of  equation  (3.2-7).  One 
of  the  more  familiar  yet  complicated  methods  is  through  the  use  of  Green's 
function  (see  Lamb  et  al  .  (1974)).  The  following  solution  is  presented 
here  because  it  is  thought  to  be  novel.  This  approach  to  the  solution 
was  arrived  at  independently,  and  was  not  discovered  during  the  extensive 
literature  search  undertaken  prior  to  writing  this  dissertation. 

The  approach  is  to  assume  a  product  solution  to  equation  (3.2-7), 
apply  boundary  condition  (3.2-7a)  to  evaluate  arbitrary  constants,  and 
finally,  verify  that  the  result  satisfies  boundary  conditions  (3.2-7b)  and 
(3.2-7c).  The  form  of  the  product  solution  is 

X(x,y,2)  =  C-F(y,x).G(z,x) 

where  C  is  an  arbitrary  constant.     Substituting  the  assumed  product  solu- 
tion into  equation   (3.2-7)  and  rearranging  yields 

Jj-T7  9E_K    ^1-      IriT^G       ,,     S^G-,  MOON 

oj  dZ 

where  it  is  understood  that  F  is  a  function  of  y  and  x  only  and  G  1s  a 
function  of  z  and  x  only.  Consider  tne  follov/ing  argument  as  applied  to 
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equation  (3.2-8).  For  a  fixed  x,  a  change  in  y,  and  hence  F,  would  in 
general  change  the  left  hand  side  of  the  equality  but  not  the  right  hand 
side.  The  opposite  is  true  for  a  fixed  x,  and  a  change  in  z,  and  hence 
G.  A  change  on  the  right  hand  side  should  occur,  but  no  change  on  the 
left  hand  side.  For  equation  (3.2-8)  to  be  a  true  equality  it  must  be 
satisfied  for  all  values  of  x,  y,  and  z  and  not  for  just  arbitrary  values. 
The  inescapable  conclusion  is  that  both  sides  of  equation  (3.2-8)  must  be 
a  constant,  such  that  coordinate  changes  leave  the  value  unchanged.  This 
constant  is  referred  to  in  mathematical  literature  as  a  separation  constant. 

It  is  expedient  at  this  point  to  consider  this  separation  constant 
zero,  rather  than  to  carry  it  through  the  manipulations  only  to  have  it 
lumped  in  with  the  other  constants  in  the  evaluation  of  the  boundary  con- 
dition (3.2-7a).  This  effectively  separates  (3.2-8)  into  the  following 
two  equations 


2  9^F  _  9F 

a 


9y 


2  dx  (3.2-9) 


and 


2 
B  -2-3^  (3.2-10) 

dZ 

2  —  2  — 

where  a    =  K  /u  and  6     =  K^/u.     Notice  that  the  two  equations  are  identi- 
cal  in  form  so  that  only  one  of  them  has  to  be  solved. 

To  solve  equation  (3.2-9),  again  assume  a  product  solution  of  the 

t 
form 


F(y,x)  =  exp(py  +  qx)  (3.2-11) 


t 
It  is  a  product  since  exp(a)exp(b)  =  exp(a+bj 
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where  exp  =  2.71828...,  the  Naperian  exponential.     Substituting  equation 
(3.2-11)  Into  (3.2-9)  yields 


or 


2  2 
a  p  F(y,x)  =  qF(y,x) 


2  2 
q  =  a  p 


2  2  2  2 

Let  p     =  -s    ,   so  that  q  =  -a  s    ,  and  p  =  -js.     Then 


2  2. 


F(y,x)  =  exp(:!"jsy  -  a  s  x) 


is  a  solution  to  equation   (3.2-9).     A  particular  value  of  s  yields  a  par- 
ticular solution.     The  total   solution  is  the  linear  combination  of  all 
possible  solutions.     The  total   solution  then,   is 


F(y,x)  =        exp(*jys  -  a  xs  )ds 


which  evaluates  as  follows: 


F(y,x)  = 


F(y,x)  = 


f  2     2 

exp(-a  xs   )[cos(ys)*j  sin(ys)]ds 


00 


2     2  r  7     9 

exp(-a  xs   )cos(ys)dsij         exp(-a  xs    )sin(ys)ds 


F(y,x)  = 


f  2     2 

2        exp(-a  xs   )cos(ys)ds 

0 


(3.2-12) 


The  last  step  follows  since  the  first  integral   is  an  even  function  while 
the  second  integral    is  an  odd  function.     From  a  table  of  comn^.on  integrals, 
equation   (3.2-12)  becomes 
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F{y,x)  =  ^  (-)    exp(-y  /4a  x) 

Since  the  partial  differential  equation  for  G(z,x)  equation 
(3.2-10),  is  identical  in  form  to  equation  (3.2-9),  we  have  for  G(z,x) 

«/   \   l/'n'\l/2    I     2, .,,2, 
G(z,x)  =  g  (-)    exp(-z  /46  x) 

Recall  that  the  form  of  the  solution  was  assumed  to  be 

X  =  C-F(y,x)-G(z,x) 

thus 

2     2 

°'   ^   ^       4a2x   4a2x 

which  upon  substitution  for  a  and  3  becomes 

= -%  (>P  [-  f^  -  ^1  (3.2-13) 

Equation  (3.2-13)  can  be  placed  in  the  standard  form  of  a  Gaussian 

function  if  we  let 
ua2 
K  =  -^  (3.2-14) 


y 


2x 


and 


K,  >-^  (3,2-14a) 


Then 


X  =  §^  exp[-  1(^)2  -  \if)h  (3.2-15) 


y  z      y     z 
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The  arbitrary  constant  C  can  now  be  determined  from  the  boundary 
condition  (3.2-7a),  which  is  repeated  here  for  convenience 


u  X  dydz  =  Q 


(3.2-7a) 


Substitution  of  equation  (3.2-15)  into  boundary  condition   (3.2-7a)  and 
rearranging  yields 


u2TrC 


^  exp[-  l(^)  ] 


f      1 


,    exp[-  h^)  ^^y    dz  =  Q 

y  y 


The  integral   in  the  brackets  evaluates  to  v^  (from  the  properties  of  a 
Gaussian  function).     This  leaves 


(2ir}  .      Cu 


^  exp[-  i(-?-)  ]dz  =  Q 


and  again,  the  integral  evaluates  to  /Zii,  and  we  have 


47r^Cu  =  Q 


or 


C  =  5Al 

2 


Substitution  of  this  constant  into  equation  (3.2-15)  yields 


""••^-i  =  2%  «pc-  |(^i'  -  hip 


(3.2-16) 
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Equation  (3.2-15)  is  the  steady  state  solution  to  the  dispersion 
equation  for  a  single  source  of  strength  Q  at  the  origin.  It  is  actually 
the  product  of  a  Gaussian  function  of  y  and  a  G=>ussian  function  of  z. 
At  first  glance  it  does  not  appear  to  be  a  function  of  x.  However, 
inspection  of  equations  (3.2-14)  and  (3.2-14a)  reveal  that  a  and  a    are 
both  functions  of  x.  The  rationale  for  substituting  the  a's  (standard 
deviations)  for  the  K's  (eddy  diffusivity  constants)  is  so  the  Gaussian 
functions  appear  in  standard  form  and  so  a  comparison  with  the  result  from 
the  statistical  approach  (section  3.3)  will  be  possible. 

For  completeness,  we  must  verify  that  boundary  conditions  (3.2-7b) 
and  (3.2-7c)  are  also  met.  This  is  easiest  to  do  by  inspecting  equation 
(3.2-13).  The  presence  of  the  exponential  makes  it  obvious  that  X 
approaches  zero  as  r  approaches  infinity,  and  so  boundary  condition 
(3.2-7b)  is  satisfied.  To  verify  that  X  approaches  infinity  as  r  approaches 
zero  (boundary  condition  (3.2-7c)),  notice  that  the  presence  of  the: 
squared  distance  in  the  exponential  numerator  means  the  exponential 
approaches  unity  as  r  approaches  zero.  Then,  the  x  in  the  denominator 
outside  the  exponential  is  the  governing  quantity,  and  it  obviously 
approaches  infinity  as  r  and  hence,  x  approaches  zero. 

Equation  (3.2-15)  is  rewritten  below  to  emphasize  the  product  of  two 
Gaussian  functions 

2  2 

X(x,y,z)  =^    — ^— exp[-  l(-f-)  ]  -J— exp[-^;f)  ]  (3.2-17) 
u  v^  a^      ^  "^z  /2Tr  a  ^   ^y 

Equation  (3.2-17)  shows  that  the  concentrations  are  proportional  to  the 
source  strength,  Q,  and  inversely  proportional  to  the  mean  wind  speed,  il 
in  the  downwind  direction,  a^  .  Also  the  concentrations  are  Gaussian  dis- 
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tributed  1n  the  vertical,  a,,  and  horizontal,  a  ,  directions. 

The  solution  obtained  above,  equation  (3.2-17),  can  be  verified  under 
more  general  conditions  than  were  used  in  the  derivation.  The  practice 
used  here  v/as  to  obtain  the  solution  by  appropriate  mathematical  manipu- 
lations without  justifying  each  steo.  In  the  end  the  solution  is  justi- 
fied on  its  own  merits;  that  is,  it  satisfies  the  partial  differential 
equation  (3,2-7),  and  boundary  conditions  (3.2-7a),  (3.2-7b),  and  (3.2-7c). 
In  the  next  section  we  present  the  statistical  approach  to  dispersion 
modeling.  We  will  see  that  it,  too,  leads  to  a  product  solution  of  two 
Gaussian  functions. 

3 . 3  Steady  State  Statistical  Model 

The  statistical  approach  to  dispersion  modeling  of  pollutants  emitted 
into  the  atmosphere  will,  as  we  will  see,  lead  to  a  Gaussian  formulation 
much  the  same  as  the  gradient  transport  theory  formulation.  We  will 
present  the  condition  under  which  the  two  lead  to  an  identical  result. 
The  statistical  approach  has  its  beginnings  in  a  study  of  turbulence  of 
the  atmosphere  by  Taylor  (1921).  Unlike  the  gradient  transport  theory, 
this  approach  does  not  begin  by  considering  the  physical  laws  of  materials 
in  the  atmosphere.  Instead,  probability  theory  is  utilized  to  evaluate 
the  likelihood  of  the  material  being  modeled  traversing  from  source  to  a 
particular  point  in  the  atmospheric  space.  The  concentration  of  the 
material  at  that  point  is  then  proportional  to  its  probability  of  arriving 
there. 

The  following  analogy  is  described  by  Slada  (1968)  and  is  included 
here  in  a  summary  version  to  emphasize  some  of  the  underlying  assumptions 
of  the  statistical  approach.  A  classroom  experiment  is  performed  by  an 
instructor  and  his  students  in  which  the  instructor  tosses  a  coin  and 
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passes  it  to  the  student  to  his  right  or  left  in  the  middle  of  the  front 
row,  depending  on  whether  the  coin  comes  up  tieads  or  tails.  That  student, 
in  turn,  continues  the  experiment  by  tossing  the  coin  and  passing  it 
over  his  right  or  left  shoulder  to  a  student  in  the  second  row,  depending 
on  the  outcome  of  the  toss.  Similarly,  the  students  continue  tossing 
the  coin  and  passing  it  back  until  it  reaches  the  back  row.  If  this 
process  is  repeated  with  a  large  number  of  coins,  a  trend  begins  to 
appear  with  most  of  the  coins  going  to  students  near  the  middle  of  the 
last  row  and  fewest  to  those  near  the  end  of  the  last  row. 

There  are  several  features  of  this  experiment  that  have  parallels  in 
statistical  atmospheric  diffusion.  First,  the  probabilistic  nature  of 
the  travel  from  source  to  final  position  (point  in  atmosphere)  as  illus- 
trated by  the  mechanism  of  tossing  the  coin.  Second,  conservation  of 
coins  (mass)  must  not  be  violated,  i.e.,  all  coins  must  be  accounted  for. 
Third,  deposition  may  occur  if  someone  puts  coins  in  his  pocket.  Fourth, 
a  distribution  of  coins  (concentrations)  is  observed  to  occur  with 
increasing  regularity  as  the  number  of  coins  injected  (emitted)  into  the 
experiment  (atmosphere)  increases. 

The  experiment  above  has  been  mathematically  formalized  by 
Chandrasekhar  (1943),  who  showed  that  the  resulting  distribution  is 
Bernoulli's  distribution.  For  large  number  of  coins  (pollutants  emitted), 
the  central  limit  theorem  requires  that  the  distribution  approach  the 
Gaussian  function.  For  "spreading"  in  the  horizontal  and  vertical 
direction  and  transport  by  the  mean  wind,  u",  in  the  downwind  direction 
(a),  this  becomes 


-X' 


-  -~—  exp[-l(^)^   ^_exp[.  l<x/j      (3.3.1) 

v"^^  ",         z     Y^Tf  a  V 
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where,  Q,  the  source  strength  and,  u,   are  as  a  result  of  the  conservation 
of  mass.  The  standard  deviations,  ct^  and  a  ,   are  statistical  properties 
necessary  to  represent  this  dispersion.  How  they  are  determined  is 
discussed  shortly.  Notice  that  equation  (3.3-1)  is  identical  to  equation 
(3.2-17).  Equation  (3.3-1)  is  only  presented  here  for  discussion  and 
comparison  with  the  solution  to  the  gradient  transport  model.  Its 
derivation  can  be  found  in  Gifford  (1955). 

Although  equations  (3.3-1)  and  (3.2-17)  appear  identical,  there  is 
a  subtle  difference  in  the  two.  Equation  (3.3-1)  is  derived  with  standard 
deviations  (a's)  which  are  statistical  in  nature.  Equation  (3.2-17)  is 
derived  with  eddy  diffusivity  constants  (see  equation  (3.2-13)).  The 
connection  between  the  two  is  through  equations  (3.2-14)  and  (3.2-14a) 
repeated  here  for  convenience 


and 


S' 

2x 

%' 

h' 

z 
2x 

(3.2-14) 


(3.2-14a) 


Many  researchers  assume  the  eddy  diffusivity  constants  are  indepen- 
dent of  position  as  their  name  implies  (see  Lamb  et  al .  (1974),  Slade 
(1968),  Tesche  et_al^  (1976)).  This  assumption  leads  to  the  following 
incorrect  result,  which  they  take  to  be  the  fundamental  relationship 
connecting  the  gradient  transport  theory  to  the  statistical  dispersion 
theory 


dt  ^^-""^  "  2K.      i  =  y  or  2  (3.3-2) 
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Equation  (3.3-2)  is  valid  only  if  the  assumption  of  the  K's  being 
independent  of  position  is  valid.  We  will  show  that  assumption  leads  to 
a  relationship  between  the  a's  and  the  x-coordinate  that  is  not  substan- 
tiated in  field  tests.  We  will  develop  the  argument  for  the  K  and  a 
relationship  (3.2-14a)  but  the  same  also  holds  true  for  K  and  a  , 

y   y 

equation  (3.2-14). 

If  K  is  not  a  function  of  x  then 

^2^  =  a^x  (3.3-3) 

2 
where  a  is  a  proportionality  constant.  This  leads  to 


-  2    2 

\f    -  u  a  X  _  a  77 


or  K^  is  proportional  to  the  mean  wind  speed,  IF.  Intuitively,  we  feel 
that  the  turbulence  in  the  z-direction  should  not  be  a  function  only  of 
the  mean  wind  speed  in  the  x-direction  since  they  are  perpendicular. 
Even  more  convincing  that  equation  (3.3-2)  is  not  true,  in  general,  is 
that  field  tests  conducted  by  Pasquill  (1961)  have  shown  that  a     follows 
a  power  law  equation  of  the  form 


0^  =  d^x  (3.3-4) 


or 


o^^  -~   aV^  (3.3-5)        f 


which  is  in  direct  opposition  to  equation  (3.3-3)  except  whenever  b=l/2. 
The  parameters  a  and  b  are  related  to  atmospheric  stability  to  be  dis- 
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ciissed  next.     Pasquill    (1962)  has  reported  values  of  b  ranging  from  0.4 
to  2.1  depending  on  the  stability  class  of  the  atmosphere. 

Assuming  equation   (3.3-5)  to  be  valid  instead  of  equation  (3.3-3), 
since  it  has  been  experimentally  verified,   then  K    must  be  a  function  of 
the  x-coordinate,  since  from  equation   (3.2-14a) 


y    _  a^u     2b-l 

S^T"^  (3.3-6) 


thus 


u 


d    ,     2,       2K    J  dK 


d    ,_  2.       .„     .    „     ^\ 


dt  (^z  )  =  ^\  +  2x  ^  (3.3.7) 

where  the  last  step  follows  since,  u"  =  dx/dt.  Equation  (3.3-7)  is  seen 
to  contain  an  extra  term  when  compared  with  equation  (3.3-2).  Substitu- 
tion of  equation  (3.3-6)  into  equation  (3.3-7)  yields 

^(,^2)^2A^2b-1^2,^A(2,_^)^2b-2 


d  /_  2, 


-2   -  2 

ua^    ua 


(a/)  =-f-  +  ^(2b-l) 


dt  ^"2  •    x     X 
2 


d  ,  2.     ^^ 
dt-  ('^z 


(a/)  =  2b 
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i  K^)  -  4bK,  (3.3-8) 

Again  we  have  two  equations,   (3.3-8)  and  (3.3-2)  that  contradict  each 
other  unless  b  =  1/2,  which  is  simply  not  experimentally  verified. 

As  mentioned  earlier,  the  gradient  transport  theory  and  the  statis- 
tical approach  lead  to  identical  results,  providing  an  appropriate  rela- 
tionship connecting  the  eddy  diffusivity  constants  and  the  standard 
deviations  is  used.     We  have  shown  that  by  choosing  the  eddy  diffusivity 
constants  as  functions  of  the  x-coordinate  (see  equation  (3.3-6)),  and 
adopting  Pasquill's   (1961)  power  law  form  for  the  standard  deviations 

^Z  =  ax''  (3.3-9) 


and 


^y  ^  ex  (3.3-9a) 

leads  to  an  appropriate  connection  in  that  the  fundamental  connection 

between  the  two  theories,  equations  (3.2-14)  and  (3.2-14a),  is  satisfied. 

We  choose  the  power  law  form  for  the  standard  deviation  terms,  a 

z 

and  a^,  since  experimental  field  tests  have  shown  it  is  an  appropriate 
form.  The  parameters  a,  b,  c,  and  d  are  functions  of  the  stability  of 
the  atmosphere.  Pasquill  (1961)  has  categorized  the  atmosphere  into  six 
classes,  A  through  F,  where  A  is  the  most  unstable  (most  turbulent)  and 
F  is  the  most  stable  (least  turbulent).  His  original  categorization  was 
based  on  wind  speed,  cloudiness,  and  sunshine  intensity.  More  recent 
models  (Roberts  et  al .  (1970),  Busse  and  Zimmerman  (1973),  and  Air 
Quality  Display  Model  (1969))  have  retained  the  same  categories,  (A~F), 
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but  use  the  vertical  temperature  gradient  as  a  means  of  selecting  the 
appropriate  stability  class.  The  stability  of  the  atmosphere  is  best 
defined  as  the  capability  for  enhancing  or  resisting  vertical  motion. 
While  the  stability  is  obviously  a  function  of  many  meteorological  parai- 
meters  (particularly  wind  gustiness  and  vertical  temperature  gradients) 
today's  models,  almost  universally,  consider  only  vertical  temperature 
gradients  (lapse  rates). 

Since  we  will  make  extensive  use  of  the  neutral  stability  class 
(D),  we  take  this  opportunity  to  describe  it.  If  a  volume  of  air  is 
forced  upward,  it  will  encounter  lower  pressures  and  expand  and  cool. 
Theoretically,  if  no  heat  is  exchanged  between  the  volume  of  air  and  the 
surrounding  atmosphere,  the  volume  of  air  will  cool  at  about  the  rate  of 
one  degree  centigrade  per  one  hundred  meters  vertical  (-TC/lOO  M),  This 
lapse  rate  is  termed  the  dry  adiabatic  lapse  rate  and  is  characteristic 
of  neutral  stability,  which  means  there  is  no  tendency  for  the  volume  of 
air  to  gain  or  lose  bouyancy. 

If  the  lapse  rate  is  greater  than  the  dry   adiabatic  lapse  rate 
(cooling  faster  than  one  degree  centigrade  per  one  hundred  meters),  the 
atmospheric  condition  is  termed  unstable  (superadiabatic) ,  and  this  class 
tends  to  enhance  vertical  motion  in  the  atmosphere  and  therefore  disperses 
more.  If  the  lapse  rate  is  less  than  the  dry  adiabatic  rate  (subadiabatic), 
the  atmosphere  is  slightly  stable,  and  a  volume  of  air  forced  upwards 
becomes  more  dense  and  will,  as  a  result,  experience  a  downward  force 
returing  it  to  its  original  height.  Stable  conditions  do  not  favor  dis- 
persion, and  this  condition  is  commonly  called  an  "inversion".  Neutral 
conditions  are  most  i-epresentative  of  long  term  or  average  conditions, 
since  they  are  associated  with  overcast  skies  and  moderate  wind  speeds 
(Turner  (1970)). 
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In  this  work  we  chose  to  use  the  gradient  transport  theory  or 
K-theory  to  derive  the  dispersion  model,  since  it  could  be  related  to 
basic  physical  laws.  It  should  be  pointed  out  that  there  are  those  who 
are  critical  of  K-theory  (see  Calder  (1965)).  A  Russian  worker,  Monin, 
(1959),  refers  to  K-theory  as  a  "semi-empirical"  theory  of  diffusion. 
Actually,  it  is,  since  Pick's  law  is  empirical  (but  then  so  is  Newton's 
law). 

The  most  compelling  reason  for  acceptance  of  the  Gaussian  formula- 
tion as  valid  is  that  a  widely  different  approach,  the  statistical 
approach,  leads  to  the  same  formulation.  The  acceptance  of  the  Gaussian 
formulation  is  so  widespread  among  researchers  that  all  the  major  dis- 
persion models  in  use  today  incorporate  the  Gaussian  distribution. 

Finally,  we  are  able  to  avoid  the  argument  over  which  of  the  two 
approaches  to  the  Gaussian  formulation  is  most  valid,  since  we  are  actual- 
ly using  a  hybrid  of  gradient  transport  theory  and  statistical  theory. 
Our  dispersion  model  is  based  on  the  solution  of  a  simplified  diffusion 
equation  and  our  model  incorporates  the  statistical  properties  of  the 
atmosphere  through  the  use  of  standard  deviations  related  to  the  x-coordi- 
nate. 

In  the  next  section  our  steady  state  dispersion  model  is  further 
modified  to  be  representative  of  long  term  (monthly  or  seasonal)  concen- 
trations. In  addition  boundary  reflections  are  considered  and  incorporated 
into  our  mathematical  dispersion  model. 

3.4  Long  Term  Dispersion  Model 

Up  to  this  point  in  our  development  of  an  appropriate  dispersion 
model,  we  have  assumed  steady  state,  i.e.,  variations  with  time  are  set 
equal  to  zero.  The  length  of  time  for  which  our  mathematical  model 
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(equation  (3.2-17))  is  valid  varies  from  a  few  minutes  to  an  hour  or  so, 
depending  on  the  length  of  time  the  mean  wind  speed  and  direction,  stability 
parameters,  and  emissions  are  constant.  One  of  the  main  contributions  of 
this  research  is  through  the  use  of  stochastic  production  costing  as 
a  basis  for  the  comparison  of  the  effects  of  one  dispatching  strategy 
relative  to  another  strategy  on  ground  level  concentrations  of  sulfur 
dioxide.  In  order  for  that  comparison  to  be  meaningful,  the  time  frames 
of  the  dispersion  model  and  the  production  costing  model  must  be  the  same. 

As  will  be  shown  in  Chapter  4,  the  appropriate  time  scale  for  the 
production  costing  model  is  a  month,  a  season,  or  even  a  year.  Thus, 
we  now  undertake  to  modify  our  steady  state  dispersion  model  to  a  long 
term  dispersion  model,  where  it  is  now  understood  that  long  term  implies 
monthly,  seasonal,  or  annual  estimates  of  the  concentrations  of  sulfur 
dioxide. 

Before  developing  the  long  term  model  for  dispersion,  it  is  con- 
venient to  modify  the  model  through  a  coordinate  translation.  Equation 
(3.2-17)  was  derived,  in  Section  3.2,  as  if  the  source  of  the  emissions 
were  at  the  origin  (0,0,0)  of  the  cartesian  coordinate  system.  Actually, 
for  electric  power  plants  (the  only  sources  considered  here),  the  emissions 
are  from  a  tall  stack  (see  Figure  1).  Inspection  of  Figure  1  reveals  that 
the  plume  actually  disperses  symmetrically  about  an  "effective"  stack 
height,  H,  rather  than  the  actual  stack  height,  H  .  The  effective  stack 
height  is  defined  as  the  height  at  which  the  centerline  of  the  plume 
becomes  essentially  horizontal  (dotted  line  in  Figure  1).  Fonnulas  for 
determining  the  effective  stack  height  are  discussed  in  Section  3.5. 
It  simplifies  matters  to  place  the  origin  of  the  coordinate  system  at  the 
base  of  the  stack.  The  source,  then,  is  at  position  (0,0,H).  Equation 
(3.2-17)  with  this  coordinate  system  translation  becomes 
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Figure  1.  Coordinate  System  Showing  the  Plume  Dispersing 
About  an  Effective  Height  in  the  Vertical  and 
Horizontal  Directions. 
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X  =  ^   -Z^-  exp  [-  -l(.|^i)2]  --L_  exp  [-  i(J^)2]       (3.4.1 ) 
u    /27TO2         "^z     /2Tra^       ^  ''y 

Throughout  this  dissertation,  we  are  most  concerned  with  ground 
level  concentrations  of  sulfur  dioxide,  and  we  take  this  opportunity  to 
consider  only  ground  level  (z=0)  concentrations.  Equation  (3.4-1)  becomes 

X  =  ^  -^  exp  [-  i(i^)2]   -1-  exp  [-  \{-f)h  (3.4-2) 

z  y       "^ 

Notice  that  equation  (3.4-2)  is  no  longer  a  function  of  z.  VJe  shall 
present  empirical,  formulas  in  Section  3.5  for  determining  H  and  IT,  so 
that  they  become  deterministic  inputs.  Thus,  the  only  variables  in 
equation  (3.4-2)  are  x,  y,  and  Q.  Remember  that  a  and  a    are  functions 
of  X.  Equation  (3.4-2)  is  rewritten  bel 


ow 


=  1 


where 


X(x,y,q)  =  ^  •  g(H)  •  g(y)  (3.4-3) 

u 


g(H)  =  -^  exp  [-  -kil-)2]  (3.4.33) 

/27ra^       "  z 


and 


g(y)  =  -=—  exp  [-  \{^)h  (3.4-3b) 

y      ^ 

Implicit  in  equation  (3.4-3)  is  that  the  x-direction  is  the  downwind 
direction.  Obviously,  for  long  periods  of  time  this  direction  will  change 
many  times.  A  particularly  convenient  way  to  represent  these  changes  in 
wind  direction  and  also  wind  speed  is  to  use  historical  wind  rose  data 
representative  of  the  time  period  for  which  concentration  estimates  are 
desired.  Usually,  15  wind  directions  are  assumed  (16  major  compass  points). 
The  wind  rose,  then,  gives  for  each  wind  direction  the  frequency  of 
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occurrence  and  the  mean  wind  speed.  For  a  good  first  order  approximation, 
the  neutral  stability  class  D  may  be  assumed  (Turner  (1970)). 

Since  16  wind  directions  have  been  assumed,  it  is  convenient  to 
think  of  the  stack  as  located  at  the  center  of  a  circle  of  radius,  x, 
and  the  circle  to  consist  of  16  sectors  with  each  sector  centered  along 
one  of  the  16  major  compass  points.  Consider,  for  example,  the  sector 
shown  in  Figure  2. 

The  wind  rose  will  give  information  about  the  frequency  of  the  wind 
and  its  speed,  for  example,  in  the  direction  of  east.  Thus,  to  obtain  a 
meaningful  estimate  of  the  concentration  as  a  function  of  x  in  the  eastern 
direction,  it  is  necessary  to  consider  the  average  concentrations  over  the 
crosswind  (horizontal)  distance  from  y=  -d/2  to  y=d/2.  We  will  define  this 
crosswind  average  concentration  later,  but  for  now  we  distinguish  it 
from  the  concentration  function  in  equation  (3.4-3)  by  placing  a  bar  over 
it.  Thus,  the  crosswind  average  concentration  function  is  xi>^,Q). 
This  crosswind  averaging  will,  in  effect,  remove  the  explicit  dependence 
of  the  concentrations  on  the  y-coordinate  and  will  make  the  concentrations 
functions  of  the  x-coordinate  (the  radial  distance  from  the  stack),  the 
emissions,  Q,  and  the  frequency,  f,  with  which  the  wind  blows  toward  that 
sector.  So,  for  a  given  value  of  x  and  a  particular  wind  rose,  15 
concentration  estimates  are  obtained  (one  for  each  direction). 

The  average  concentration  estimate  for  a  particular  direction  as  a 
function  of  the  radial  distance,  x,  is  the  product  of  the  crosswind 
average  concentration,  J,   and  the  percentage  of  time  the  wind  blows  in 
that  direction  (the  frequency,  f,  in  that  direction).  Thus 

x(x,Q,f)  -  f  -  x(x,Q)  (3.4-4) 
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Figure  2.  Eastern  Sector  Showing  Downwind  and 
Crosswind  Directions. 
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It  is  important  to  remember  that  the  concentration  function,  x, 
(equation  (3.4-3))  is  a  function  of  only  two  random  variables,  Q  and  y, 
for  a  given  value  of  x.  All  other  parameters  are  deterministic.  Then, 
by  definition 


x(x,Q)  =     x(x,y,Q)p(Q,y)  dQdy 


(3.4-5) 


where  p(Q,y)  is  the  joint  probability  density  function  of  Q  and  y. 
For  obvious  physical  reasons,  the  emissions,  Q,  are  independent  (statis- 
tically) of  the  y-coordinate  and  so 


p(Q>.y)  =  p(Q)p(y) 


and,  then 


X(X:,Q) 


g(H) 


Qg(y)p(Q)p(y)dQdy 


y^ 


(3.4-6) 


which  can  be  written  as 


x(x,Q)  =  ^  1  g(y)p(y)  f  Qp(Q)dO  dy 

U    J  Y  JQ 


(3.4-7) 


Note  that  the  expression  in  the  brackets  in  equation  (3.4-7)  is,  by 
definition,  the  average  or  mean  value  of  emissions.  We  shall  label  it  Q. 
Then,  equation  (3.4-7)  becomes 


X(x)  =  ^g(H)    g(y)p(y)dy 
u     •'y 


(3.4-8) 


where  the  functional  dependence  on  Q  has  been  dropped,  since  Q  is 
deterministic.  How  Q  is  computed  is  discussed  in  Section  3.5. 

It  is  reasonable  to  assume  that  a  concentration  at  a  particular 
distance,  y,  from  the  x  axis  within  the  sector  is  as  likely  as  a 
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concentration  at  any  other  value  of  y  within  the  sector  for  the  same  value 
of  the  x-coordinate.  That  is,  within  the  sector  and  a  particular  radial 
distance,  x,  from  the  stack,  the  concentrations  should  be  weighted  equally. 
Then 

p(y)  =  c 


where  c  is  a  constant.  But  application  of  the  definition  of  a  probability 
density  function  requires 


d 
2 

d 
2 


c  dy  =  1 


or 


c  =  1/d 


This  yields  for  x(x) 


-(  \     Q  g  H) 

u     ^-  ^  

2    y 

Consider  the  change  of  variables 


2  _J r  ■I/y^2-,^ 

d  -^r-  exp  [-  j{-^)   ]dy 

/2TTa..         ^y 


S=  y/o^ 


a. 


then  dS  =  ^^,  since  a,  is  not  a  function  of  y.  So 

y      J 


x(x)=?^f'^-L 


d   1     _  exp  (-  i  s^)dS 


where 


Si  =  d/2ay 
From  Figure  2,  d  is  given  by 

d  =  xe 
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(In  fact,  d=2  tan(T^)x  =  .398x,  whereas  d-  -^y.   -  .393x  or  about  a  one 

lb  10 


percent  error) .     So 


^M'^^m-Gis^) 


where  g 

1      '  ^■' 


/2tF 
and 


1     r2, 


G(S^)  -  -^  exp  (-  ^  S'')dS 


■h 


n     _       2ttx     _       X 


1       16-2a         16a 

y        y 

For  the  assumed  neutral   stability  class  D  (Pasquill    (1961)) 

a^  Vcx^  =  0.32  x°-^S 

so 

S^   =  0.614  x°-^2 

Refer  to  Table  2  for  several  values  of  x,  S-,,   and  G(S,)  evaluated 
from  standard  mathematical  tables.  For  distances  from  the  stack  of  one 
kilometer  or  more,  the  function  G(S, )  is  for  all  practical  purposes  equal 
to  unity.  Then 

x(x)  =  2^±g(H) 

u 
or 

{2-nr^   u  >^^2       2  a^ 
Finally,  from  equation  (3.4-4),  we  have 

X(x.f)   =  — 16        _fSL  g^p  [-_  l(iL)2j  (3_4_9j 

(2.)^'^^  uxa^  2a^ 
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Table  2.  Integrated  Gaussian  Function  for  Several 
Values  of  Downwind  Distance. 


Downwind  Limit  of  Integrated 

Distance  Integration,  Gaussian 

in  Meters,  S.  Function, 

X '  G(S^) 


100  1.69  .925 

500  2.41  .982 

1000  2.81  .994 

5000  4.00  .999 

>5000  >4.00  1.000 
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Equation  (3.4-9)  is  yery   similar  to  one  developed  by  Calder  (1971), 
who  used  jDolar  rather  than  rectangular  coordinates. 

Equation  (3.4-9),  when  modified  to  include  boundary  reflections, 
will  be  the  long  term  dispersion  model  used  throughout  the  remainder  of 
this  dissertation.  In  Chapter  5  the  results  of  Chapters  3  and  4  are 
applied  to  the  study  of  a  realistic  size  electric  power  system  with 
multiple  generating  units.  Also,  a  computer  program  developed  during  this 
research  is  presented  and  discussed.  For  distances  closer  to  the  source 
than  one  kilometer,  equation  (3.4-9)  must  be  modified.  This  is  done  also 
in  Chapter  5 . 

Equation  (3.4-9)  is  plotted  in  Figure  3  as  a  function  of  the  radial 
distance,  x.  It  is  seen  to  begin  at  zero  concentration,  increase  to  a 
maximum,  and  then  fall  off  rapidly  to  zero  again. 

As  mentioned  in  Chapter  2,  the  EPA  has  set  a  radius  of  liability 
of  25  miles  (40.2  km)  for  each  point  source.  It  is  instructive  to 
determine  from  equation  (3.4-9),  at  what  radial  distance,  x,  the  maximum 
concentration  occurs  as  a  function  of  effective  stack  height  for  a  given 
set  of  meteorological  conditions  and  emissions.  Equation  (3.4-9)  is 
rewritten  below  as 

X=k  (--^)exp[-  1(^)2] 

A 

where  k  is  a  constant  containing  the  fixed  emissions  and  meteorological 
parameters.  If  x  is  differentiated  with  respect  to  x  and  the  result  set 
equal  to  zero,  the  maximum  is  found  to  occur  at 


1 


1 
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5     10     15    20     25 
Downwind  Distance  in  Kilometers 


Figure  3.  Plot  of  Concentrations  Vs.  Downwind 
Distance  from  the  Source. 
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Figure  4  shows  a  plot  of  equation  (3.4-10)  as  a  function  of  H  from  one 
hundred  meters  to  one  kilometer.  The  values  of  a  and  b  are  0.22  and 
0.78  respectively  (Pasquill  (1961)).  Figure  4  reveals  that  the  maximum 
occurs  farther  downwind  for  increasingly  higher  effective  stack  heights. 
Even  for  the  unreal istically  high  one  kilometer  effective  stack  height, 
the  maximum  occurs  well  within  the  25-mile  liability  radius  set  by  the  EPA. 
One  must  keep  in  mind  that  while  higher  effective  stack  heights  produce 
maximums  farther  downwind,  these  maximums  are  decreasing  significantly  with 
higher  effective  stacks.  This  is  shown  in  Figure  5. 

Our  long  term  dispersion  model,  equation  (3.4-9),  is  valid  only  if 
the  effective  stack  height  is  far  enough  away  from  the  earth's  surface 
that  its  presence  can  be  ignored.  Physically  this  can  never  be  the  case, 
and  we  must  include  the  effect  of  the  earth's  surface  in  the  model.  The 
most  effective  way  to  include  this  physical  barrier  is  to  assume  total 
reflection  of  the  pollutant  at  the  earth's  surface.  While  some  deposition 
at  the  ground  is  known  to  occur,  we  are  being  conservative  in  our  estimates 
by  ignoring  it.  The  most  expedient  method  for  handling  these  reflections, 
mathematically,  is  through  the  concept  of  a  virtual  source  or  image  source 
located  symmetrically  with  respect  to  the  ground  plane,  to  the  actual 
source.  This  technique  is  widely  used  in  heat-conduction  theory. 

The  translation  of  the  source  from  origin  (0,0,0)  to  effective  stack 
height  (0,03H)  caused  the  vertical  dispersion  term  to  change  from  z  to  z-H. 
For  a  (virtual)  source  located  a  distance  underground  of  H  meters,  the 
vertical  dispersion  due  to  this  second  source  should  disperse  about  -H, 
or  appear  in  the  exponential  as  a  z+H  term.  Thus,  the  vertical  dispersion 
term  becomes 

(vertical  dispersion)  =  — ^—  exp  [-  h^h  +  exp   [-  i(^^] 

/27Ta        ^  "^z  "^    "^z 
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Figure  4.  Plot  of  Downwind  Distance  to  Maximum  Concentration 
vs.  Effective  Stack  Height  of  Source. 
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Figure  5.  Plot  of  Maximum  Concentration  vs.  Effective 
Stack  Height. 


■"■^'^l^'-rf ' f^fWf ^^flW^it-i"*— T^«*- ■T**—ril*-^y^'WlW^si.-H:^" 'IwMfatf  f  — ws^?  "f  y  -^^m'tr^ik 


(-JlM«r-^r '  iJM-il ilK£.=^=ic<- *.-,^r->  r 


MitmS^v^-rjstisnim 


66 


Reflections  of  the  pollutant  may  also  occLir  off  a  stable  layer  of 

air  above  the  stack  (inversion)  in  much  the  same  way.  If  the  height 

of  the  inversion  (sometimes  called  mixing  height)  is  H  ,  then  we  can 

include  its  effect  by  assuming  a  virtual  source  at  !-i  +  (H  -H)  =  2  H  -H, 

mm       m 

Then  the  vertical  dispersion  exponents  have  the  following  three  terms: 

z-H  (due  to  the  actual  source),  z+H  (due  to  ground  level  reflections), 

and  Z-2H  +H  (due  to  inversion  layer  reflections).  There  can  also  be 

reflections  off  the  inversion  layer  after  reflections  off  the  ground.  This 

leads  to  a  term,  z-2H  -H.  Likewise,  there  can  be  reflections  off  the 

ground  after  reflections  off  the  inversion  layer.  This  leads  to  a  term, 

z+2Hj^-H.  Consider  the  case  where  reflections  occur,  first  off  the  ground, 

second  off  the  inversion,  and  third  off  the  ground.  This  leads  to  a 

term,  z+21-1  +H.  The  symmetry  is  complete  now,  since  all  possible  combinations 

of  signs  in  front  of  H  and  H  have  been  shov/n.  Additional  multiple 

reflections  are  theoretically  possible  and  each  new  reflection  changes 

only  the  numerical  coefficient  on  H  .  This  is  because  of  the  additional 

m 

distance  traversed  from  ground  to  inversion,  H  .  Consideration  of  all 

m  


possible  reflection  combinations  leads  to  the  following  infinite  set  of 
terms 


z-H,  I     [z-(2nH  +H)],  I     Cz-(2nH  -H)] 
n=l       ^  n=l       ^ 


z+H,  I     [z+(2nH  +H)],  I     [z+(2nH  -H)] 
n=l      ^  n=l      ^ 


It  is  important  to  remember  that  these  terms  actually  appear  in  exponentials, 
and  we  have  a  sum  of  exponentials  rather  than  the  sum  of  the  terms  above. 

Since  we  are  primarily  concerned  with  ground  level  concentrations, 
we  let  z=0  and  after  some  grouping  of  like  terms,  we  have 
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Equation  (3.4-s3a)  is  now  changed  to  include  reflections  and  is  written 

g(H)  =— ^gR(H)  (3.4-11) 


Our  long  term  dispersion  model  including  reflections  is  written 
below  in  its  entirety  for  future  reference 


^^;p7?V^   exp[-2(5^)]. 

»        ,  2nH^-H  ,        ,  2nH  +H  , 
I       exp[-  i(^L)2]  .  exp[-  i(-JL_)2]        (3.4-12) 
n- 1  z  z 

In  the  computer  program  discussed  in  Chapter  6,  it  is  seldom 
necessary  to  include  terms  beyond  n=3,  since  the  large  exponents  cause 
the  exponentials  to  approach  zero  very   rapidly.  An  internal  check  in 
the  program  compares  values  of  the  infinite  summation  terms  to  the  first 
term  (n=0  term)  beginning  with  n=l .  Values  that  do  not  contribute  more 
than  0.1  percent  are  discarded. 

Inversion  layers  aloft  are  not  typical  of  long  term  meteorological 
conditions.  The  most  useful  long  term  dispersion  model,  then,  includes 
ground  level  reflections  but  not  inversion  layer  reflections.  This  can 
be  accomplished  in  equation  (3.4-12)  by  setting  H^^  =  »,  which  causes  the 
infinite  summation  to  vanish,  leaving 

v^Y  f\  -  2.03fQ    r  "!/  H  ^2. 

^'  '  ^  '  "^  ^^p!^"  ?V^  -^  (3.4-13) 

Z  "  z 

Equation  (3.4-13)  is  the  one  we  will  use  most  frequently  in  modeling 
the  pollutant  dispersion  from  electric  power  plants. 
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In  the  next  section  the  plume  rise  formulas  are  presented.  A  formula 
for  calculating  u   is  presented,  and  the  techniques  for  calculating  gross 
emissions  is  developed. 

3.5  Plume  Rise,  Wind  Speed  at  Stack  Height,  and  Emissions  Formulas 
We  begin  this  section  with  a  discussion  of  the  many  plume  rise 
formulas  in  use  today.  Inspection  of  Figure  1  shows  that  there  are  two 
distinct  regions  of  the  plume  dispersion.  In  the  first  region,  the  plume 
rises  vertically  from  the  stack,  due  to  its  stack  exit  velocity,  to  some 
higher  position  termed  the  effective  stack  height.  From  there  the 
prevailing  mean  wind  carries  the  plume  downwind  while  it  disperses.  It 
is  important  to  remember  that  even  though  the  plume  enters  the  atmosphere 
at  (0,0,H^),  the  actual  stack  height,  it  is  treated  mathematically  as  if 
it  entered  the  atmosphere  at  (0,0,H),  the  effective  stack  height  with  zero 
vertical  velocity. 

The  early  attempts  to  estimate  the  effective  stack  of  a  buoyant  plume 
were  primarily  empirical  formulas  based  on  observation.  The  point  at 
which  the  plume  centerline  becomes  essentially  horizontal  is  highly 
subjective  and  explains  why  there  have  been  so  many  different  plume  rise 
formulas  developed  over  the  years.  Tesche  et  al .  (1976)  quotes  Briggs 
(1969)  as  saying,  "There  were  over  30  documented  formulas  as  of  1969, 
and  he  [Briggs]  estimated  that  two  more  would  be  added  each  year," 

Disappointed  by  the  vast  differences  (sometimes  as  much  as  three  to 
one)  in  estimates  of  plume  rise,  researchers  began  attempting  to  develop 
plume  rise  models  using  rigorous  physical  laws.  Two  general  approaches 
developed  and  it  was  encouraging,  at  first,  that  both  led  to  the  same 
"two- thirds  power  law."  The  first  approach,  proposed  by  Morton  et  al , 
(1956),  uses  entrainment  theory,  -rfhile  the  second  approach.  Csanady  (1955), 
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uses  gradient-transfer  closure  arguments.  We  will  outline  the  entrainment 
theory,  pointing  out  the  assumptions  that  affect  its  application  to  our 
dispersion  model . 

The  approach  begins  by  equating  the  buoyant  force  of  the  plume  to  the 
time  rate  of  change  of  the  vertical  momentum  of  the  plume 

Trr^gp  Y  "  ^  (Tir^pw)  (3.5-1) 

where 

r  =  radius  of  plume  (m) 

2 

g  =  acceleration  of  gravity  (m/sec  ) 

3 
p  =  density  of  plume  material  (ug/m  ) 

T  =  ambient  temperature  (°K) 

AT  =  T  -T  =  temperature  excess  between  plume  and 

surroundings  (°K) 

w  =  vertical  plume  speed  (m/sec) 

The  plume  is  assumed  to  be  adiabatic  which  leads  to 

2— 

Qu  =  irr  uAT  =  constant  (3.5-2) 

where 

3 

Q  =  heat  emission  (m  °K/sec) 
n 

u  =  mean  wind  speed  throughout  tha  plume  in  the  downwind 

direction  (m/sec) 

For  equation  (3.5-2)  to  be  strictly  valid,  requires  that 

AT  a  — p 
T 

since  u  is  assumed  constant. 
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Using  equation   (3.5-2),  we  can  replace  the  left-hand  side  of 
equation  (3.5-1)  with  pgQu/Tu.     The  right-hand  side  can  be  rewritten  as 

dz 
Recognizing  that  ^  =  w,  we  have 

d    /     2     \  d    /   2   \ 

w  -^  (rrr  pw)   =  'fipw  ^  (r  w) 

which  becomes 

<i    f     2     .  2  dw  ^        2  dr" 

w  ^  (.irr  pw)  =  TTpwr    -^  +  TTpw    ■^— - 

Equation  (3.5-1)  can  now  be  rewritten  as 

^,2^  dw  ^   2  dri  ^  9^1  , 

dz      dz   TlT  ^"^-^  "*' 

Briggs  (1971)  assumes  that  the  plume  radius  increases  linearly  with 
height  during  the  ascent  to  effective  stack  height,  i.e. 

r  =  Y2 

where  y  is  a  constant  determined  experimentally.  Now  equation  (3.5-3) 
becomes 


2_2  .  dw   ,       2  2,_       ^ 

Tu 


irYZw^+7TYw2z  =  -^jizr  (3.5-4) 


Assume  w=az  ,  where  a  and  n  are  constants,  as  the  solution  to 
equation  (3.5-4).  This  yields 

2  2  n   n-1   „  2  2  2n    ^% 

■vry  z   az  naz    +  2tp,-  a  z  z  =  —tr 

Tu 

or 

,,2  2  2n+l  ,  „  2  2  2n+l   ^^H  ,,  ,.  ,. 

na  Try  z    -^  2a  Try  z    =  -~  (3.b-5) 

Tu 
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Because  of  the  adiabatic  assumption,  the  right-hand  side  of  equation  (3.5-5) 
is  a  constant  (independent  of  z).  Thus,  we  must  have 

2n+l  =  0 


or 


Then 


or 


n.-i 


1  2  2  ^  „  2  2   9Ql| 

■5-  a  Try  +  2a  tty  =  — z: 
^  Tu 


Stty  Tu 


so 


2gQ    1/2  _,,„ 

Stty  Tu 


where  A  is  recognized  as  a  constant. 

Since  by  definition,  the  downwind  speed  is 

—   dx  , 

•^  =  dt 

and  the  vertical  wind  speed  of  the  plume  is 

'^  "  dt 
their  ratio  is 

iz  =  w  =  w^zl  (3.5-7) 

u    u 

since  by  equation  (3.5-6)  w  is  a  function  of  z.  Rearranging  equation 
(3.5-7)  and  integrating,  we  have 
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^        u 


A  3/2--- 


or 

z=  (3AX)2/3^  (|a2)1/3(^)2/3 

2ir     ^      u 
Substituting  equation  (3«5-6)  and  grouping  terms,  we  have 

3gQn   1/3  ^  ,,^ 
z  =  [-^]    [^]"'^  (3.5-8) 

2Try  Tu      u 

Equation  (3.5-8)  shows  that  the  height,  z,  the  plume  rises  above  the 
actual  stack  height  is  proportional  to  the  downwind  distance,  x,  raised 
to  the  two-thirds  power;  hence,  it  is  often  called  Briggs'  two-thirds  power 
law. 

Equation  (3.5-8)  is  often  written  in  terms  of  the  buoyancy  flux,  F, 
which  is  defined  as 

F  =  gvr^  ^-  (3.5-9) 

where  v  is  the  velocity  of  the  plume  material,  and  the  other  terms  are  as 
previously  defined.  Note  for  the  initial  buoyancy  flux,  v=w  (the  stack 
exit  velocity),  while  for  the  buoyancy  flux  at  or  near  the  effective  stack 
height,  v=u  (the  mean  downwind  speed). 

Using  equation  (3.5-9)  we  can  rewrite  equation  (3.5-8)  as 

2y  v     u 

It  should  be  pointed  out  that  while  equation  (3.5-10)  is  valid  in  form  for 
all  x,  it  is  of  use  to  us  only  near  the  stack  (v=w)  or  near  the  effective 
stack  height  (v=u),  since  v  is  unknown  elsewhere.  This  is  net  a  limitation. 
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actually,  since  for  our  dispersion  model  we  are  only  interested  in  z  at 
or  near  the  effective  stack  height. 

Baker  and  Jacobs  (1971)  have  suggested  that  the  plume  has  essentially 
reached  its  effective  height  at  a  downwind  distance  of  10  times  the  actual 
stack  height.  With  this  modification  and  v^u",  equation  (3.5-10)  becomes 

,.  (.3^)1/3  i  „    2/3  (3^^_„j 


o  -  a 

Zy  u 


where 


c     2-  AT 
F  =  gr  u  ^ 


Remembering  that  the  adiabatic  assumption  required  constant  heat  emission, 

Q^,  and  therefore  M  a  -^  ,   we  note  that  the  buoyancy  flux,  F,  is  also 

t         "^ 
a  constant.   Since  it  is  easier  in  practice  to  determine  the  parameters  at 

stack  height,  we  chose  the  initial  buoyancy  flux  rather  than  the  final. 

This  means 

F  =  gr  w  4j.  (3.5-12) 

where  it  is  understood  that  r  is  now  the  inside  radius  of  the  actual 
stack,  and  w  is  the  stack  gas  exit  velocity.  Tesche  et  al .  (1976)  report 
that  the  entrainment  constant,  y,   is  approximately  0.66  for  a  neutral 
atmosphere. 

To  avoid  confusion  with  the  z-coordinate,  in  the  future  we  will  define 
the  increase  above  the  actual  stack  as  Ah.  Thus  equations  (3.5-11)  and 
(3.5-12)  become 

Ah'^z.i!^(M:^)V3,„,2/3  ,  ^ 

u 


'This  is  equivalent  to  assuming  conservation  of  momentum. 
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and  the  effective  stack  height  is,  then,  given  by 

H  =  Hg  +  Ah  (3.5-14) 

While  there  are  literally  dozens  of  plume  rise  formulas  in  the 
literature,  we  have  chosen  equation  (3.5-13)  for  our  model,  since  it  is 
partly  based  on  physical  laws.  In  addition,  a  second  theoretical  approach, 
the  gradient- transfer  closure  (Csanady  (1965)),  leads  to  the  same  form. 
It  should  be  remembered  that  the  model  is  empirical  due  to  the  entrainment 
constant  relationship  (r  =  yz) .     The  assumption  of  adiabatic  conditions  is 
also  a  limitation,  although  not  a  serious  one  for  neutral  atmospheric 
conditions  (Liu  (1975a)). 

In  the  interest  of  being  objective,  we  should  point  out  that 
Thomas  et  al .  (1970)  have  studied  several  plume  rise  models  in  field 
tests  including  the  one  presented  here,  and  report  that  it  has  a  slight 
tendency  to  underpredict  the  rise  for  light  winds  and  overpredict  for 
strong  winds.  They  report,  however,  that  "the  use  of  the  2/3  power  law 
formula  is  considered  preferable,  provided  information  for  the  meteorological 
parameters  is  available."  We  again  emphasize  that  since  our  decisions  will 
be  based  on  comparisons,  a  good  consistent  estimate  of  plume  rise  is 
satisfactory  for  our  purposes. 

Recently,  there  has  evolved  a  new  approach  to  modeling  the  plume  rise 
phenomenon.  This  approach  involves  numerical  solutions  to  equation  (3.1-6), 
the  generalized  dispersion  equation,  near  the  stack,  (see  Shir  (1970) 
and  Liu  (1975b)).  This  approach  requires  extensive  computer  computations 
and  to  date  has  not  been  shown  to  be  a  significant  improvement  over  the 
two-thirds  law  model  we  chose.  This  approach  shows  much  promise,  however, 
and  may  in  the  future  replace  all  the  empirical  models.  It  is  considered 
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too  early  in  its  developmental  stage  to  be  of  use  to  us  now,  but 
undoubtedly,  this  will  change  in  the  future. 

The  mean  wind  speed,  u,  is  used  in  both  the  long  term  dispersion 
model  and  the  plume  rise  model.  Wind  speed  is  known  to  increase  with 
height,  i.e.,  u=u(z).  Ideally,  the  mean  wind  speed  for  the  dispersion 
model  should  be  representative  of  the  wind  about  the  plume  center! ine 
(effective  stack  height),  while  the  mean  wind  speed  for  the  plume  rise 
model  should  be  representative  of  the  wind  throughout  the  vertical  region 
from  the  actual  stack  height  to  the  effective  stack  height.  To  avoid  this 
complication,  most  models,  including  the  one  used  here,  use  the  mean  wind 
speed  at  stack  height,  H^.  This  allows  for  sampling  the  wind  at  one  height 
and  eliminates  the  errors  that  would  arise  from  inaccurate  determinations 
of  the  effective  stack  height. 

For  the  plume  rise  model,  the  entrainment  coefficient,  y,   is 
determined  empirically  from  field  tests  where  u"  is  selected  at  the  actual 
stack  height.  Likewise,  for  the  dispersion  model  the  stability  parameters, 
a  and  b,  are  determined  empirically  with  u"  selected  at  the  actual  stack 
height.  Since  u  at  10  meters  off  the  ground  is  the  most  common  method  for 
reporting  mean  wind  speeds,  we  need  a  formula  to  relate  the  mean  wind 
speed  at  stack  height  to  the  mean  wind  speed  at  10  meters. 

Busse  and  Zimmerman  (1973)  suggest  a  power  law  relation  for  the 

neutral  atmosphere  of  the  form 

H 
"^=  "^lO  ^TU)^'''^  (3.5-15) 

where  u-jq  is  the  mean  wind  speed  at  10  meters,  and  the  other  parameters 
are  as  previously  defined. 
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Finally,  we  complete  this  section  by  presenting  the  formula  for 
determining  gross  emissions  of  sulfur  dioxide  from  the  stack.  Several 

definitions  will  be  presented  first. 

t 
Heat-input  can  be  defined  as  the  rate  of  heat,  in  BTU's/hr  necessary 

to  produce  a  given  amount  of  electrical  power  out  of  the  generator. 

Several  of  the  emissions  standards  set  by  the  U.  S.  Government  are  in 

terms  of  pounds  per  million  BTU   heat-input  (see  Chapter  1).  In 


mathematical  terms,  the  heat- input,  Hj,  is 


H,  =  P  ■  H^ 


where 


P  =  electrical  power  generated  in  megawatts,  (MW) 
H  =  heat  rate  in  BTU's  per  10^  watt- hour  (BTU/KWH) 

As  an  example,  a  500  MW  output  for  a  unit  having  a  9,000  BTU/ 
heat  rate  requires 

Hj  =  ( 500x1 0^w)( 9000  |^) 

Hj  =  4.5x10^  BTU/H 

or  four  and  one-half  billion  BTU's  per  hour. 

If  this  were  a  coal-fired  unit,  the  coal  consumption,  C  ,  would  be 

"l 
where  H  is  the  heating  value  of  the  coal  in  BTU's  per  pound.  For  this 


We  regret  this  sudden  shift  from  the  Mi(S  system  of  units  to  the  British 
system  of  units,  but  the  use  of  BTU's  is  so  widespread  in  the  electric 
power  industry  that  the  use  of  Joules  or  Watt-seconds  would  be  distracting. 

+4" 

Even  the  U.S.  Government  is  inconsistent,  since  it  uses  the  MKS  system 

3 

for  concentrations  (ug/i^  )■ 
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unit  using  a  coal   having  12,000  BTU's  per  pound,  the  coal   consumption 
would  be 

P     „  4.5x10^  BTU/H 
c  "     12000  BTU/lb 

C^  =  3.75x10^  Ib/H 

or  375,000  pounds  per  hour  of  coal. 

If  this  coal  contained  p  percent  sulfur  then  the  sulfur  con- 
sumption S  would  be 

C  •  p 
S  =  "■      ^ 


'c    100 

where  p^  is  the  sulfur  content  of  the  coal  in  percent.  For  this  unit, 
if  p  =  three  percent,  then  the  sulfur  consumption  is 

•3! 

(3.75x10^  Ib/H)  •  3 
^C  "        100 

S  =  1.125x10^  Ib/H 
c 

or  11,250  pounds  per  hour  of  sulfur. 

During  the  combustion  process,  the  sulfur  is  oxidized  to  sulfur 
dioxide  (SO2).  The  atomic  weight  of  sulfur  is  16  while  the  two  oxygen 
atoms  together  have  an  atomic  weight  of  16.  Thus  the  weight  of  sulfur 
dioxide  is  twice  that  of  elemental  sulfur.  This  means  that  for  each  pound 
of  sulfur  consumed  in  the  furnace,  there  will  be  two  pounds  of  sulfur 
dioxide  produced  and  emitted  from  the  stack.  This  assumes  complete 
(100  percent)  oxidation.   Thus,  for  this  unit,  the  sulfur  dioxide 
produced,  Q,  is  given  by 


t 
The  EPA  allows  certain  furnaces  95  percent  oxidation,  or  a  weight  factor 

of  1.9  rather  than  2.0.  Our  choice  of  2.0  is  in  keeping  with  a  con- 
servative estimate. 
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Q  -  2  ■  S^ 

Q  =  2  •  1.125x10^  Ib/H 

Q  =  2.250x10^  Ib/H 

or  22,500  pounds  per  hour  of  sulfur  dioxide  produced. 

It  is  desirable  to  have  a  formula  for  calculating  gross  emissions 

of  sulfur  dioxide,  Q,  directly.  This  is  shown  below: 

20P  •  H  •  p 
Q  =  _ s  (3.5-16) 

where  all  of  the  parameters  are  as  defined  above.  Substitution  of  P=500, 

H^=9,000,  p^=3,  and  H^=l 2,000,  verifies  the  result  as  22,500  pounds  per 

hour  of  sulfur  dioxide  produced. 

Since  the  U.S.  Government  has  selected  micrograms  per  cubic  meter 

for  reporting  concentration,  equation  (3.5-16)  is  modified  to  the 

following  form 

(2.53xlO^)P  •  H  •  p 
Q  =  H,   '    '  -     (3.5-17) 

where  Q  is  now  given  by  micrograms  per  second.  As  is  shown  in  Chapter  6, 
it  is  necessary  to  have  both  equations  (3.5-16)  and  (3.5-17).  Ground 
level  concentrations  are  determined  through  the  dispersion  model  using 
equation  (3.5-17),  after  it  is  modified  in  Chapter  4  to  reflect  average 
(long  term)  emissions.  Emission  standards,  on  the  other  hand,  are  computed 
on  the  basis  of  pounds  of  sulfur  dioxide  per  million  BTU  heat-input, 
which  requires  the  use  of  equation  (3.5-16). 
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CHAPTER  4 
THE  METHOD  OF  PRODUCTION  COSTING 

In  this  chapter  we  present  the  method  of  production  costing  as  it 
is  used  by  today's  electric  utility  companies  for  estimating  the  long 
term  fuel  costs  of  the  supply  of  electrical  energy.  We  will  extend  the 
theory  to  include  other  functions  of  the  output  of  the  generating  units. 
In  particular,  we  will  be  interested  in  the  application  of  this  extension 
to  the  estimation  of  long  term  emissions  of  sulfur  dioxide. 

The  production  costing  method,  then,  is  applied  to  the  problem  of 
determining  the  best  estimates  of  long  term  emissions  of  sulfur  dioxide 
from  each  source  (generating  unit),  and  these  estimates  are  used  as 
inputs  to  the  mathematical  dispersion  model  presented  in  Chapter  3.  Through 
the  use  of  the  production  costing  method  and  the  dispersion  model, 
decisions  can  be  made  based  on  the  order  of  dispatching  generating  units 
that  minimize  the  incremental  deterioration  of  ground  level  air  quality 
throughout  some  appropriate  geographical  region. 

In  the  past  the  electric  utilities  have  experienced  exponential 
growth  rates  due  to  the  ever- increasing  demand  for  electric  energy.  A 
doubling  rate  of  ten  years  has  not  been  uncomm.on,  and  some  areas,  for 
example  Florida,  have  experienced  doubling  rates  closer  to  seven  years. 
Faced  with  these  short  lead  times,  utility  system  planners  have  relied 
heavily  on  larger  and  larger  capacity  generators.  Along  with  these 
increased  capacities  has  naturally  been  an  increased  complexity  in  design 
and  operation  of  these  large  units.  As  manufacturers  of  electric  generators 
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pushed  the  limits  of  their  technologies,  a  new  phenomenon  began  to 
emerge--that  of  an  increasing  unavailability  of  these  units  due  to  forced 
outages  brought  about  by  malfunctions  in  boilers  or  auxiliary  parts. 
It  became  obvious  during  the  early  sixties  that  a  technique  for  incor- 
porating these  forced  outage  rates  into  the  various  system  planning  models 
had  to  be  developed. 

In  1967,  French  researchers  (Baleriaux  et  al .  (1967))  introduced  a 
probabilistic  simulation  technique,  that  included  forced  outage  rates 
for  estimating  the  fuel  costs  of  operating  various  electric  generating 
units  within  an  electric  utility  system.  Hydroelectric  and  pumped-storage 
projects  were  assigned  fictitious  fuel  costs  based  on  operating  costs  for 
purposes  of  comparison.  Booth  (1972)  incorporated  a  dynamic  programming 
algorithm  with  the  method  introduced  by  Baleriaux  et  al .  (1957)  and 
introduced  his  generation  expansion  program  as  applied  to  the  State 
Electricity  Commission  of  Victoria,  Australia. 

In  addition  to  consideration  of  the  forced  outage  rates  of  the 
generating  units,  the  production  costing  model  as  applied  by  Booth  (1972) 
could  consider  all  types  of  generating  units,  i.e.,  nuclear,  fossil, 
hydroelectric,  gas  turbine,  combined  cycle,  and  pumped-storage.  This  has 
gained  the  production  costing  method  increasing  acceptance  within  the 
electric  utility  industry  in  several  areas  of  application.  Sager  et  al . 
(1972)  have  also  used  the  method  for  evaluating  generation  production  costs, 
Sullivan  (1974)  has  used  the  method  for  generation  reserve  planning, 
while  Sullivan  and  Hilson  (1975)  have  used  it  in  estimating  ground  level 
concentrations  of  sulfur  dioxide  from  nearby  electric  generating  plants. 

In  Section  4.1  we  will  extend  the  mathematical  theory  of  production 
costing  to  include  any  monotonically  increasing  or  decreasing  function  of 
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the  output,  P,  of  a  generator.  Since  in  Chapter  3  we  were  able  to  v/rite 
the  emission  rate  of  sulfur  dioxide  (equation  {3.5-17))  as  a  function  of 
the  output,  P,  we  can  make  use  of  this  extension  of  the  theory  to  predict 
long  term  emission  rates  from  individual  generators.  In  order  to  derive 
this  extension,  it  will  be  necessary  to  take  a  detailed  look  at  the 
mathematics  of  production  costing. 

4.1  The  Load  Duration  Curve 

At  the  heart  of  production  costing  are  three  important  concepts: 

1.  Load  duration  curves 

2.  Forced  outage  rates 

3.  Order  of  dispatching  generating  units 

In  this  section  we  will  discuss  the  use  of  a  load  duration  curve  in 
estimating  the  long  term  output  of  a  generator.  We  will  derive  the 
necessary  mathematics  to  extend  its  use  to  any  monotonically  increasing 
or  decreasing  function  of  the  output. 

One  of  the  major  contributions  of  this  research  is  the  extension  of 
the  production  costing  method  to  include  these  generalized  functions  of 
the  output  of  individual  generators.  The  only  restriction  is  that  these 
functions  be  monotonically  increasing  or  decreasing  functions  of  the 
output,  P.  Simply  stated,  this  means  there  must  be  a  one  to  one  relation- 
ship between  the  function  of  P  and  the  output,  P.  Thus,  for  a  given  value 
of  P  there  must  be  a  unique  value  of  the  function  of  P  and  vice  versa. 
Emission  rates,  fuel  costs,  and  fuel  consumption  are  functions  of  the 
output,  P,  that  satisfy  this  restriction.  In  practice,  then,  this  is  a 
restriction  only  in  the  mathematical  sense. 

Figure  6  shows  a  typical  load  duration  curve.  Load  duration  curves 
are  normally  assembled  front  historical  data.  The  abscissa  represents 
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T/2 


Tirrje  in  Hours 


Figure  6.     Typical   Load  Duration  Curve. 
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the  number  of  hours  that  the  system  load,  P.  ,  equaled'''  or  exceeded  the 
ordinate  value  of  load.  Thus,  if  T  were  the  total  time  period  under 
consideration,  then  for  T/2  hours  the  system  load  equaled  or  exceeded 
P    p  megawatts  (see  Figure  6).  If  the  values  along  the  abscissa  are 
divided  by  the  time  period,  T,  the  abscissa  values  become  a  fraction  of 
time.  In  this  case  we  can  say  that  one  half  the  time  the  system  load 
equaled  or  exceeded  P  _,„  megawatts,  or  the  system  load  equaled  or 

L  ,  1  /  ^ 

exceeded  P,  -p,„  with  a  probability  of  0.5.  Note  that  by  considering  the 
abscissa  values  as  probabilities,  we  can  use  historical  data  as  repre- 
sentative of  future  data  for  the  purpose  of  predicting  future  loads,  pro- 
viding that  the  load  duration  curve  for  the  time  period,  T,  is  a  valid 
representation  of  the  future  period  (see  Jenkins  and  Joy  (1974)). 

Generally,  the  shape  of  the  load  duration  curves  for  different  periods 
of  time  are  the  same.  The  actual  values  of  the  loads  may  be  different. 
Most  electric  utilities  have  a  fair  degree  of  expertise  in  load  forecasting, 
and  in  this  work,  we  assume  that  the  load  duration  curves  are  representative 
of  the  particular  future  time  period  under  consideration. 

In  the  development  of  probabilistic  simulation  models,  it  has  been 
more  convenient  to  reverse  the  roles  of  the  ordinate  and  abscissa  of  load 
duration  curves.  Figure  7  shows  this  "inverted"  load  duration  curve. 
The  interpretation,  now,  is  that  the  ordinate  value  is  the  forecasted 
probability  that  the  abscissa  value  (load)  is  equaled  or  exceeded  during 
the  time  period  for  which  the  curve  is  valid.  From  this  point  on  when  we 
speak  of  a  load  duration  curve,  it  should  be  understood  that  we  mean  an 
inverted  load  duration  curve. 


''"Theoretically,  one  should  only  include  the  number  of  hours  Pl  exceeded  the 
ordinate  value  of  load.  In  practice  it  is  convenient  to  include  also  the 
number  of  hours  Pl  equaled  the  ordinate  value  of  load.  The  distinction  is 
important  only  for  discontinuous  load  duration  curves  which  are  not  used 
here  (Sullivan  (1975)). 
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Figure  7.     Typical    "Inverted''  Load  Duration  Curve, 
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Since  we  will  be  deriving  the  long  tern  estimate  of  functions  of 
the  generator  output,  P,  from  fundamental  mathematics,  it  is  necessary 
to  develop  the  load  duration  curve  from  more  basic  probability  theory. 

Figure  8  shows  a  typical  probability  density  function  of  system 
loads,  fi_(Pi_)  >  where  the  "L"  indicates  system  load.  The  cumulative 
probability  distribution  function,  F'i_^^L^'  ^^   defined  as  (see  Figure  9a) 


F'l(Pl)  ^ 


fP. 


^\i\U\ 


The  backwards  cumulative  distribution  function,  F  (P  ),  can  be  defined 
as  (see  Figure  9b) 


Equation  (4.1-1)  can  be  rewritten  as 

rP, 


(4.1-1) 


Fl(P,) 


f,  (P,  )dP  -   ^f,  (P,  )dP, 


(4.1-2) 


where  the  first  integral  on  the  right  hand  side  of  equation  (4.1-2) 
is,  by  definition  of  the  probability  density  function,  equal  to  one. 
Then,  equation  (4.1-2)  may  be  written  as 


f.iPO 


fL(PL)dPL 


(4.1-3) 


Inspection  of  Figure  9b  reveals  that  it  is  identical  to  the  (inverted) 
load  duration  curve  of  Figure  7.  Equation  (4.1-3),  then,  is  the 
mathematical  description  of  a  load  duration  curve. 


In  Chapter  3  we  used  p(-)  for  the  probability  density  function.  From 

now  on  we  use  f(-)  for  the  probability  density  function  to  avoid  confusion 

with  the  generator  outputs,  P,  and  the  system  load,  P,  . 
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System  Load  in  Megawatts 


Figure  8.  Typical  Probability  Density  Function  of 
System  Loads. 
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In  order  to  illustrate  the  technique  for  determining  the  expected 

output  of  a  generator  from  the  load  duration  curve,  consider  the  following 

argument.  Refer  to  Figure  10.  Suppose  that  the  system  generators  are 

scheduled  such  that  generator  1  is  loaded  until  its  maximum  before 

generator  2  is  loaded,  and  then  generator  3  follows  after  generator  2 

is  completely  loaded,  and  so  forth. 

Suppose  also  that  the  maximum  output  of  generator  1  is  P, 

1  ,max 

which  is  less  than  P.   ■  ,  the  system  minimum  forecasted  load.  Then 
generating  unit  1  must  be  loaded  to  its  maximum  for  the  entire  time  period, 
since  from  the  load  duration  curve  (Figure  10),  F,  (P,  )=1.0,  or  the 
probability  that  the  system  load  exceeds  the  capability  of  unit  1  is 
equal  to  unity. 

In  this  case  the  output  of  unit  1  is  deterministic  rather  than 
random.  The  probability  density  function  of  a  deterministic  value  is  a 
dirac-delta  function.  Thus 

where  6(-)  is  the  notation  for  the  dirac-delta  function.  Then  the 


expected  value  of  the  output  of  unit  1  is  by  definition 

P 

_  A  r ""  (-1  ,max 

P  =   P,f(PjdP  =  I    P^flPjdP  (4.1-5) 

where  the  last  step  follows  as  a  practical  consideration.  Substitution 

or  equation  (4.1-4)  into  (4.1-5)  and  evaluating  yields 
P, 


^1  = 


1  ,max 

P,  6(Pt-Pt  )dPT 
„  1  M  l,max'  1 
u 


or 


P  =  P 
1    l.max 
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Figure  10.  Load  Duration  CiirvG 
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where  the  last  step  follows  from  the  sampling  theorem  for  dirac-delta 

functions. 

It  is  no  surprise  that  the  expected  value  of  unit  1  turned  out  to 

be  its  maximum  capability.  The  fact  that  unit  1  was  dispatched  first, 

and  its  maximum  capability  was  less  than  the  minimum  system  load, 

determined  that  it  be  fully  loaded  for  the  entire  period.  Notice  that 

the  area  under  the  load  duration  curve  in  Figure  10  from  0  to  P,  ^^„ 

I  jmax 

is  P,    ,  since  the  ordinate  value  is  unity. 
1  ,max  -^ 

Carrying  this  illustration  one  step  further,  suppose  that  it  has 

been  determined  that  generating  unit  i  will  be  loaded  only  after  the 

first  i-1  units  have  been  dispatched.  How  this  is  actually  determined  will 

be  discussed  in  Chapter  5.  This  is  the  same  as  saying  that  unit  i  will 

not  be  loaded  until  the  system  load  reaches  Pj  megawatts,  where  Py  is 

the  sum  of  the  outputs  of  the  first  i-1  units.  In  addition,  as  the  system 

load  grows  from  P-r  to  Py+P,-  „.„>  unit  i  will  be  loaded  from  0  to  P.  ^,„ 
^         111 jmax  1 ,max 

(see  Figure  11 ) . 

In  effect,  the  probability  of  dispatching  unit  i  is  F.  (P,),  the 

probability  that  the  system  load  equals  or  exceeds  P,.  Likewise,  the 

probability  that  unit  i  is  fully  loaded  is  F,  (P_  +  P.    ).  Similarly, 
^  ^  ■J  L  T    1  ,max 

the  probability  of  unit  i  being  called  on  to  supply  P  megawatts 

(0  <  P,  <  p.    )  is  the  probability  that  the  system  load  if  P^+P  or 
a    1 5 max  i  a 

greater,  which  is  precisely  F,  (P^+P  ),  Since  P  is  an  intermediate 

L  i  a        a 

output  between  no  load  and  full  load  for  unit  i,  we  know  the  probability 

of  any  output  for  unit  i  between  no  load  and  full  load.  The  average  value 

of  the  output  of  unit  i,  then,  could  be  determined  by  multiplying  each 

possible  output  by  the  probability  it  is  required  and  summing  all  these 

products.  In  the  limit,  this  average  value  for  the  output  of  unit  i  is 

simply  the  area  under  the  load  duration  curve  frorii  P..  to  P^+P. 
'^  ■'  I     T  1  ,max 
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Figure  11.  Load  Duration  Curve  Showing  Commitment 
of  Unit  i . 
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The  important  point  in  the  above  discussion  is  that  by  deciding 
beforehand  not  to  allow  unit  i  to  be  loaded  until  the  system  load  grew 
beyond  P^,  and  to  load  unit  i  to  its  maximum  as  the  system  load  grew  to 
'^T'*"'^i  max'  we  have  forced  unit  i  to  assume  the  backward  cumulative 
distribution  function  or  load  duration  curve  of  the  system  loads. 

In  the  above  discussion,  it  was  shov>/n  in  an  heuristic  manner  that 
the  expected  value  of  unit  i  was  given  by  the  area  under  the  load  duration 
curve  where  unit  i  was  placed  (see  Figure  11).  This  will  be  proven  in 
the  following  discussion. 

In  order  to  make  the  discussion  more  general,  consider  that  it  is 
required  to  find  the  expected  value  of  some  function  (not  necessarily 
linear)  of  the  output  of  a  particular  generating  unit.  The  only  restriction 
on  the  function  will  be  that  it  be  a  monotonically  increasing  or  decreasing 
function  of  the  output,  as  discussed  earlier.  In  addition,  suppose  it 
is  necessary  to  know  the  expected  value  of  the  function  over  an  interval 
(segment)  from  a  to  b.  Let  the  function  be  given  by 

Q  =  Q(P)"^ 

where  P  is  the  output  of  the  generating  unit. 

By  definition,  the  expected  value  of  this  function  over  the  interval 
a  to  b  is 

-   f^ 

Q  =    Q  f(Q)dQ  .  (4.1-6) 

■'a 
where  f(Q)  is  the  probability  density  function  of  Q.  Since  Q  is  a 


The  derivation  that  follows  is  for  some  particular  generating  unit,  e.g., 
the  i-th  unit.  However,  to  keep  the  notation  from  becoming  too  cumbersome, 
we  will  not  use  the  subscript  i,  but  it  should  be  understood. 


fii^MtV -HpHtraan^itfwm^iffH >■  w >ii tarfi»wi^BM-itiffca5Ca'a;aF»  itfin-  ii  wr-nr-i  ' 


93 


function  of  P,  equation  (4.1-6)  can  be  rewritten  as 
rb 


Q(P)f(P)dP  (4.1-7) 


where  f(P)  is  the  probability  density  function  of  the  output,  P,  of 
the  generator.  By  definition,  f(P)  can  also  be  v^ritten  as 


t-fp)  -  dPiPl  _  ::dFiP) 

'\^l  rl[)         HP 


dP      dp 

where  F'(P)  is  the  cumulative  distribution  function  and  F(P)  is  the 
backward  cumulative  distribution  function  of  P  (see  equation  (4.1-1)). 
Equation  (4.1-7)  can,  then,  be  written  as 
fb 


Q  =  - 


Q(P)dF(P)  (4.1-8) 

a 


Note  that  in  the  formulation  of  equation  (4.1-8),  the  expected  value  of 
Q  is  a  Riemann-Stijlets  integral,  since  it  is  integration  with  respect  to 
a  function  F(P),  instead  of  a  variable.  Now  suppose  that  the  unit  is 
forced  to  follow  the  load  duration  curve  from  P_  to  P-.+b-a  megawatts. 
That  is,  the  unit  is  dispatched  from  a  to  b  megawatts,  only  as  the  system 
load  grows  from  P-j.  to  P-j.+b-a  megawatts.  Then  for  the  interval, 
Pj  -  Pj  -  Pj+b-a,  for  the  system  load,  there  is  a  corresponding  interval, 
a  -  P  ^  b,  for  the  generating  unit.  These  corresponding  limits  require 
that 

P  =  Pj_-P^+a  (4.1-9) 

and 

dP  =  dP^ 

Now  since  the  generating  unit  follows  the  load  duration  curves  from 
P-j-  to  P-j-+b-a,  its  backward  cumulative  distribution  function  F(P),  may 
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be  expressed  in  terms  of  the  load  duration  curve  as 


F(P)  =  [y(PL-Pj)  -  u(Pj_-P-^-b+a)]F|_(P^) 


where  y(-)  is  the  step  function  and  is  related  to  the  dirac-delta 
function  by 


(4.1-10) 


dP 


y(P-c)  =  6(P-c) 


where  c  is  a  constant.  Then  by  the  equation  above,  we  have  for 
equation  (4.1-10) 

dF(P)  -  [5(PL-Py)-6(PL-PT-^^^)^  FL(PL)dPL 


+  [li(PL-PT)-V'(PL-PT-t>+a)]  dF^(P^) 


(4.1-11) 


Substitution  of  equation  (4.1-11)  into  (4.1-8)  and  changing  the 


limits  to  reflect  the  change  of  variables  (equation  (4.1-9) )yields 
Q  =  -!      Q(PL-PT+a)  [6(P^_-P^)-6(PL-Py-b+a)]  F|_(PL)dP^_ 


P,+b-a 


P,+b-a 


Q(P^-P.^+a)[u(P^-P^)-u(PL-P^-b+a)]  dFjP^) 


T 


The  first  integral  can  be  simplified  by  the  sampling  theorem,  and  the 

second  integral  can  be  simplified  by  noting  that  the  step  functions  are 

identical  to  the  limits  of  integration.  Hence,  they  are  unnecessary. 
This  leaves 


Q  =  -Q(a)F^(P^)+Q(b)F^(P.,.+b-a)  - 


P  +b-a 


Q(PL-P^+a)dF^_(P^)   (4.1-12) 


Then  by  use  of  the  partial  integration  formula,  the  integral  in  equation 
(4.1-12)  becomes 
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P,+b-a 
r  I 


Q(P^_-P^+a)dF^(P^)  =  -Q(b)FjP^+b-a)+Q(a)Fj_(P^) 


P-[.+b-a 
+  f    F^(PL)dQ(P^-P^+a) 


(4.1-13) 


Substitution  of  equation  (4.1-13)  into  (4.1-12)  and  simplifying  yields 


Q  = 


P^+b-a 


FL(PL)dQ(PL-PT+^) 


(4.1-14) 


The  restriction  that  Q  and  P  have  a  one  to  one  relationship  satisfies  the 
requirement  of  the  partial  integration  formula  that  Q  be  of  bounded 
variation  on  the  interval  a  to  b. 

Assume  for  the  moment  that  the  function  Q  is  the  output  P  of  the 
generator.  That  is 

Q(P)  =  P 
so 


Then 


Q  =  P 


dQ.  =  dP 


and  by  equation  (4.1-9) 
dP  =  dP 

and  so  equation  (4.1-14)  becomes 

Q  =  P 


P^+b-a 
r  T 


Fj_(PL)tiPL 


which  by  definition  is  the  area  under  the  load  duration  curve  from  Pj 
to  P-p+b-a,  as  earlier  hypothesized. 


t-jiy>c;ic:*.s;=; 
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In  Chapter  3  the  emission  rate  of  sulfur  dioxide  (in  yg/sec)  was 
shown  to  be  (equation  (3.5-17)) 


Q  = 


(Z-SSxIO^^)?  •  H^  •  p^ 


For  realistic  size  generating  units,  the  heat  rate,  H  ,  is  a  function  of 
the  generator  output,  P.  Thus  the  emission  rate  should  be  written  as 


Q  =  k  •  P  •  H^(P) 


(4.1-15) 


where  k  is  independent  of  P  and  is  given  by 
(2.53xl0^)p^ 


H. 


Differentiating  equation  (4.1-15)  with  respect  to  P  yields 


f  ^k^[P-H,(P)] 


(4.1-16) 


The  derivative  on  the  right  side  of  equation  (4.1-15)  is,  by  definition, 
the  incremental  heat  rate,  IHR(-)»  of  the  generating  unit  (in  BTU/KWH). 
Returning  to  equation  (4.1-14)  and  using  equation  (4.1-9),  we  have 


P.|.+b-a 


F^(Pj_)dQ(P) 


which  by  equation  (4.1-16)  becomes 
Q 


P^+b-a 


F^(P|_)k  IHR(?)dP 


Since  P=P|^-Py+a,  and  dP=dP|  ,  we  finally  have 

P,.+b-a 


Q  =  k 


P^ 


IHR(P|-P^+a)!-^(F^)dP, 


(4.1-17) 
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Equation  (4.1-17)  is  the  expected  value  cf  the  emission  rate  of  sulfur 
dioxide  from  a  generator  dispatched  over  the  interval  a  to  b,  corresponding 
to  the  system's  load  duration  curve  from  Pj   to  P^+b-a.  Note  that  for 
the  limits  of  integration  in  equation  (4.1-17)  the  incremental  heat  rate 
of  the  generator  assumes  values  from  a  to  b  as  would  be  expected,  k'hen 
Q  is  multiplied  by  the  time  period,  T,  (in  seconds)  for  which  the  load 
duration  curve  is  valid,  the  result  is  the  expected  emissions  of  sulfur 
dioxide  (in  micrograms)  the  generating  unit  would  produce  over  the  interval 
a  to  b.  This  concept  is  used  in  Chapter  5  to  compare  the  effect  various 
generating  units  have  on  ground  level  concentrations  of  sulfur  dioxide, 
as  each  additional  segment  of  energy  under  the  load  duration  curve  is 
dispatched. 

The  general  approach  taken  in  deriving  equation  (4.1-17)  makes  it 
easy  to  compute  the  expected  value  of  other  functions  of  P.  For  example, 
suppose  one  desires  the  expected  cost  per  hour  of  operating  unit  over 
some  range  (a  to  b)  of  its  output,  P.  Let  Q(P)  be  replaced  by  C(P),  the 
cost  per  hour.  Next  suppose  this  unit  is  placed  under  the  load  duration 
curve  from  P^  to  P,+b-a,  then  from  equation  (4.1-14) 


C  = 


P  +b-a 

F^_(PjdC(PL-Py+a) 


The  cost  per  hour  of  operating  the  unit,  considering  only  the  fuel  costs, 
is  given  by 

C(P)  =  fj,  •  P  •  H^(P)  (4.1-18) 

where  f  is  the  fuel  cost  in  dollars  per  1000  BTU's  (S/lO'^  BTU)  and  the 
other  parameters  are  as  previously  defined. 


ai>>^    fTlMBl    II   i""    i'rT«ii«'iT--'i  "   i'"" 
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Differentiation  of  equation  (4.1-18)  yields 


^  ^  -c  ■  ^  [P-H^(P)]  =  f^-IHR(P)  (4.1-19) 


Again,  using  the  change  of  variables,  P=P, -Pj+a,  we  have 
C  =  f 


P  +b-a 

IHR(P^-P^+a)F|_(P^)dP^  (4.1-20) 

Note  equations  (4.1-20)  and  (4.1-17)  differ  only  by  a  constant.  We  can 
relate  them  as 

Q  =  (k/f^)  C  (4.1-21) 

In  Chapter  5  we  use  equation  (4.1-20)  to  compute  the  costs  of  dispatching 

the  units  over  various  segments,  and  then  use  equation  (4.1-21)  to  compute 

the  emissions.  Note  for  a  nuclear  unit  or  a  hydroelectric  unit  k=0 

(since  p.=0)  and  the  emissions  are  zero.  For  this  reason,  we  cannot 
s 

compute  all  the  costs  from  the  emissions,  which  is  why  we  did  not  write 
equation  (4,1-21)  as 

C  =  (f^/k)  Q 

At  the  beginning  of  this  section  we  listed  three  important  concepts 
that  are  fundamental  to  the  production  costing  method.  The  first  item, 
the  load  duration  curve,  has  been  discussed  in  this  section.  The  second 
item,  the  forced  outage  rates,  is  considered  in  the  next  section.  The 
final  item,  the  order  of  dispatching  the  generating  units,  is  taken  up 
in  Chapter  5. 

All  of  the  results  of  this  section  are  valid  only  if  the  generating 
unit  is  available  all  of  the  time,  i.e.,  a  forced  outage  rate  of  zero. 


Obviously  this  can  never  be  the  case,  and  our  results  must  be  modified 
to  include  non-zero  forced  outage  rates.  This  is  done  in  the  next 
section. 

4 . 2  The  Effective  Load  Duration  Curve 

In  this  section  we  need  to  extend  the  results  of  the  previous  section 
to  include  non-zero  forced  outage  rates  of  the  generating  units.  The 
follovn'ng  discussion  is  similar  to  the  one  given  by  Jenkins  and  Joy  (1974), 
except  they  considered  constant  heat  rate  curves.  We  will  not  make  this 
restriction  and  so  our  results  will  be  more  general. 

Refer  to  Figure  12,  where  the  units  have  been  placed  under  the  load 

duration  curve  in  the  numerical  sequence,  1,2,3,4 Equation  (4.1-20) 

gives  the  cost  per  hour  of  operating  a  unit  over  an  interval  a  to  b,  under 

the  load  duration  curve  from  P^  to  P_+b-a.  Initially,  we  will  be 

interested  in  the  first  unit  to  be  dispatched,  which  from  Figure  12  is 

unit  1.  We  are  also  interested,  in  this  section,  only  in  the  complete 

loading  of  a  unit.  In  the  next  section  we  again  consider  segmental 

dispatching.  In  this  case,  a  equals  zero  and  b  equals  P,    ,  the 

1  ,max 

4- 

maximum  output  of  unit  1.  The  cost  per  hour  of  unit  1  is  then  given  by 
P  +P 

Ci  =fci  jp   ihRi(Pl-Pt)^(Pl)^Pl  (^-2-1) 

T 

where  Py  is  the  total  amount  of  pov/er  dispatched  ahead  of  unit  1,  which 
in  this  case  is  zero. 

As  mentioned  before,  equation  (4.2-1)  is  valid  only  for  the  period 
of  time  that  unit  1  is  available  for  supplying  power.  If  unit  1  is 


'In  Section  4.1  we  were  able  to  derive  the  results  without  the  necessity  of 
unit  number  subscripts.  In  this  section,  however,  since  we  are  considering 
several  units  we  must  use  unit  number  subscripts,  e.g.,  P-  is  the  output 
of  the  i-th  unit. 
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Figure  12.  Load  Duration  Curve  Showing 
Unit  Commitments. 
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unavailable  (forced  out  of  service)  then  all  the  remaining  units  in 
Figure  12  must  slide  to  the  left  under  the  load  duration  curve  as  shown 
in  Figure  13.  During  this  time  unit  1  cannot  generate  any  power,  since 
it  is  out  of  service.  Let  p.  be  the  availability  of  unit  i  and  q.  be 
the  unavailability  (forced  outage  rate)  of  unit  i.  Then  since  the  unit 
can  only  be  in  one  of  two  situations  (states)  we  have 

p.  +  q.  =  1 

where  p.  and  q.  are  expressed  as  fractions  less  than  one. 

The  expected  cost  of  unit  1  is  less  than  the  value  in  equation 
(4.2-1)  by  a  factor  p, ,  or 


^1  =  Pi  ^cl 


Pt+Pi ,max 

IHR^(P^_-P^)F^(P^)dPL  (4.2-2) 

^T 


Notice  that  since  units  2,  3,  and  4,  etc.,  are  dispatched  after  unit  1, 
whether  they  are  available  or  not  does  not  affect  the  position  of  unit  1, 
and  thus  does  not  affect  the  cost  given  in  equation  (4.2-2).  ■ 

Inspection  of  Figures  12  and  13  reveals  a  different  situation  for 
unit  2.  Its  cost  is  seen  to  depend  not  only  on  its  availability  but 
also  on  the  availability  of  unit  1.  This  is  because  during  the  fraction 
of  the  time  period,  T,  that  unit  1  is  unavailable,  q, ,  unit  2  is  shifted 
to  the  left  and  hence  occupies  a  different  position  under  the  load 
duration  curve. 

The  cost  for  unit  2,  then,  will  consist  of  two  terms.  The  first 
term  will  depend  on  the  availability  of  unit  1,  while  the  second  term 
will  depend  on  the  unavailability  of  unit  1.  Both  terms  depend  on  the 
availability  of  unit  2,  since  its  cost  is  zero  while  it  is  unavailable. 
The  expression  for  the  first  term  (Figure  12)  becomes 


<rf-rii-|  -~n— ■^iTTiTiriil-wil    iir  ¥   r  r  TIi  -i  ri  i  ■■  aatMafci' i        -'T?;.'  *-'■.'.-"','  — .-ir.'^—  -- .i..-^-.  .-.^-t-^.-r..-.-  ,-».,.^.^-.-^ 
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Figure  13.  Load  Duration  Curve  Showing  Unit 

Commitments  if  Unit  1  is  Unavailable. 


•stSHTTk  ^-«*atai 


••  fr7gi'fni';»'iif1'T^Vi'fn'iiYii'TTrr^-||-|  '=*r].^  »^m 


lo: 


Cg  (first  term)  -  P]P2^"c2 


V^2,niax 


(4.2-3) 


where  P  is  now  the  total  amount  of  power  dispatched  ahead  of  unit  2 

and  is  P,    .  Notice  that  equation  (4.2-3)  is  similar  in  form  to 
1  ,max 

equation  (4.2-2),  with  the  additional  parameter,  p, . 

For  the  second  term,  in  determining  the  expected  cost  of  unit  2; 
unit  1  must  be  unavailable.  In  this  case  the  cost  is  given  by 
(see  Figure  13) 


C„  (second  term)  =  qip2f  o 


2, max 
IHR2(PL)F^(P^)dP^ 


0 


(4.2-4) 


We  can  rewrite  equation  (4.2-4)  by  changing  the  limits  of  integration, 
such  that 


C„  (second  term)  =  qip2f  o 


P    +P 
1 ,max  2, max 

IHR-(P-P,    )F,  (P  -P,    )dP,     (4.2-5) 


1  ,max 


2  L  1 ,max  L  L  1  ,max   L 


Since  Pj  is  the  total  amount  of  power  dispatched  ahead  of  the  unit  under 
consideration,  and  we  are  now  considering  unit  2,  P-j.  must  equal  P-, 
Then  equation  (4.2-5)  becomes 


.max 


T  2, max 
^2  (^^^°"^  term)  -  q^p^f^^   |   IHR2(Pl-Pt)^(Pl-' 1  ^max^^^L 


(4.2-5) 


where,  in  the  interest  of  conformity  with  later  results,  we  have  left  the 

load  duration  term  in  equation  (4.2-6)  in  terms  of  P, 

^  1  ,max 

The  total  expected  cost  of  unit  2  is  the  sum  of  the  first  and 


second  terms,  and  becomes 
T  2, max 


C,   •  ,,f,,   J   IHR^IP^-P^)  [»lfL<PL^*^lF|.(\-Pl,r„x>^='fL 


{'■■2-7) 
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Comparing  equation  (4.2-2)  and  (4.2-7)  reveals  that  the  term  in  the 
brackets  of  equation  (4.2-7)  has  replaced  the  load  duration  curve  term 
in  equation  (4.2-2).  For  this  reason  it  is  often  called  the  "effective" 
load  duration  curve.  We  can  define  the  effective  load  duration  curve  as 

EFl(Pl)|^=P,F,(PlH,,F^(P,-P,_„J  (4.2-8) 

where  EF, (P, )   is  the  effective  load  duration  curve  after  one  unit  (the 

*-  *-  1 
first  unit)has  been  dispatched.  Substitution  of  equation  (4.2-8) 

into  (4.2-7)  yields 

T  2, max 
^2  "  P2^c2  I  IHR2(Pl-Pt)EF|_(P^_),   dP^  (4.2-9) 

pj  n 

where  it  should  be  remembered  that  P,  is  the  total  amount  of  power 

dispatched  before  unit  2.  In  this  case  P^  equals  P-,    .  Since  Pj 

includes  the  total  power  dispatched  ahead  of  unit  2,  whether  available 

or  not,  it  is  often  termed  the  effective  or  equivalent  load.  Simply 

stated,  the  equivalent  load  is  the  actual  load  plus  the  additional  load 

on  the  remaining  units  due  to  the  unavailability  of  units  already 

dispatched  (often  called  capacity  on  outage). 

Our  formulation  is  unique,  in  that  it  is  unnecessary  to  make  this 

arbitrary  definition  of  the  new  random  variable,  capacity  on  outage. 

We  are  able  to  avoid  this  additional  complication  through  the  change  of 

limits  of  integration  from  equation  (4.2-4)  to  (4.2-5).  Figure  14  shows 

a  plot  of  F,  (P,  ),  F,  (P,  -P.  ^_),  and  EF,  (P,  ),  .  Notice  that  the  effective 
L  L    L  L   I  ,max        L  >-  1 

load  duration  curve  makes  the  system  load  look  larger.  This  is  a  direct 

result  of  the  inclusion  of  the  forced  outage  rates,  q.,  of  the  generating 

units. 
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Figure  14.  Load  Duration  Curve  [F, (P, )],  Shifted 

Load  Duration  Cui-ve  [F,  (P, -Pn    )],  and 
Effective  Load  Duration  Curve  [EF  (P.  )]  . 
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We  can  show  this  mathematicany  by  the  following  argument.  Let 

P,    =  A,  then  we  want  to  show 
1  ,max 

p^F,  (P|_)+q^F^(P^-A)  ^  F^_(P^)  (4.2-10) 

The  following  steps  are  self-explanatory. 
q^FjP^-A)  ^  (l-p-|)Fj_(P^) 
F^P^-A)  ->  F|_(P^)  (4.2-11) 

Inequality  (4.2-11)  is  always  true, and  therefore  so  is  inequality 
(4.2-10),  since  the  backward  cumulative  distribution,  F.  {•),   never 
increases.  Thus,  a  shift  in  load  backwards  by  A  can  never  decrease. 
Notice  that  in  Figure  14,  the  shift  to  the  left  in  load  by  A  (equal  to 
Pi  m=v)  T"  effect  shifts  the  F,  (P,  )  curve  to  the  right  as  is  shown  by 

t  )mdX  L   L 

''■(P, -Pi   )• 

L  L  1 .max 

Inspection  of  Figure  13  reveals  that  the  expected  cost  of  unit  2 
is  independent  of  whether  or  not  the  units  to  the  right  (units  3,4,...) 
are  available.  Unit  3,  however,  will  depend  on  the  availability  of  both 
units  1  and  2,  since  the  forced  outage  of  either  one  or  both  shifts 
unit  3  to  the  left  under  the  load  duration  curve.  Thus  the  cost  of 
unit  3  will  have  four  terms  representing  the  four  possible  combinations 
of  units  1  and  2  being  available  or  not  available  for  producing  energy. 
These  ternis  are  similar  to  those  in  determining  the  expected  cost  of  unit  2. 
When  both  units  1  and  2  are  available,  C^  vn"ll  depend  on 


PiP2Fl(Pl) 
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If  unit  1  is  unavailable,  C  will  depend  on 

qi¥L^VPl,max) 
and  if  unit  2  is  unavailable  but  unit  1  is  available,  C"-,  will  depend  on 

Finally,  if  both  units  1  and  2  are  unavailable,  Co  will  depend  on 

^l¥L(V^l,max-''2,max) 

Similar  to  before,  we  will  define  the  effective  load  duration  curve 
after  two  units  are  dispatched  as 


E\(PL)|2-PlP2^L(PL)^qiP2^(PL-Pl,max) 

^  Pi¥l(  VP2,max)^^1^2f^L^PL-Pl  ,.ax-P2,.ax) 


Then  the  expected  cost  for  unit  3  becomes 


T  3, max 


C,  =  p.f 


3  c3 


IHR3(Pl-P^)EFl(P^)   dP, 


T 


(4.2-12) 


where  P-j-  is  the  total  amount  of  power  dispatched  ahead  of  unit  3  (to  the 

left  of  unit  3  under  the  load  duration  curve),  and  in  this  case  is  equal 

to  P,    +P^ 
1 ,max  2, max 

Using  arguments  similar  to  those  used  in  developing  the  expected 

costs  for  units  1,  2,  and  3,  we  can  write  the  expected  cost  for  the  i-th 

unit  to  be  loaded  as 


P^+P. 
T  1 ,max 


i  ^  Pi^ci  Jp  ihr,(Pl-Pt)ef,_(Pl) 


dP, 


(4.2-13) 


^.=iiiUa/!i,si3Si££Siu.t  .'^^•t^ 
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where  Pj  is  the  sum  of  all  the  power  already  dispatched,  or  is 

i-1 

T  -i^     J. max 

and  EF,  (P,  ).    is  the  effective  load  duration  curve  after  the  dispatching 
"-  "-  ji~l 

of  i-1  units,  and  is  given  by 

EF,  (P,  )    =  (P.P.P^  ...  P,  JF  (P,  ) 

LLj_-|       i     l.    O  1-1    Li_ 

*  ^'i^'>2''3  ■■■   l'i-l'\<''L-''l,n,ax' 

+   ... 

+  {^^^^^^   ...  qi.i)FL(PL-Pi,, ax- --Pi-Umax)   ^^-^-^^^ 

where  the  p's  and  q's  in  the  parentheses  above  are  the  2^  ~   possible 
combinations  of  the  i-1  units  being  available  or  unavailable. 

The  effective  load  duration  curve,  equation  (4.2-14),  is  often 
termed  a  convolution  equation.  This  is  because  equation  (4.2-14)  can  be 
derived  from  basic  statistical  mathematics,  involving  the  integral  of 
the  product  of  two  density  functions,  which  is  a  convolution  equation. 
Then  one  speaks  of  "convolving  in"  the  units  with  the  load  duration  curve 
and  "deconvolving  out"  the  units  when  they  are  removed  (unavailable). 
We  have  not  taken  this  approach,  since  it  requires  detailed  mathematics 
that  tend  to  obscure  the  results.  In  the  approach  we  have  taken,  it  is 
never  necessary  to  deconvolve  any  units.  Our  approach  includes,  at  each 
step,  the  appropriate  units  and  thus  avoids  deconvolution  with  its  well- 
recognized  problems. 


-^lll-iM)u-M*^h    in" ~r  TTTi 1"— ' ■  n  —       T^'V  -,'ii'^7-'7 ■-.    r      ^.u  .•!.  •— ^r  '-     ii^t- i^matmi ^Mfj^isixai ^a^y Wiii '    ar^ »»ul~ '  leniat^-sa — =jia»^;=^wr'c.-iBai«Mii^r^illwf=ay^^.apTi<;~'>tfij 
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In  the  next  section,  we  will  extend  the  results  of  this  section  to 
include  dispatching  segments  of  units  rather  than  the  entire  unit.  This 
makes  the  approach  amenable  to  the  comparison  process  to  be  discussed 
in  Chapter  5. 

4.3  The  Effective  Load  Duration  Curve  for  Segments  of  Units 

In  the  previous  section,  the  forced  outage  rate  of  the  i-th  unit 
was  given  by  q.,  and  was  described  as  the  fraction  of  the  time  period,  T, 
that  the  i-th  unit  was  expected  to  be  unavailable  for  supplying  energy. 
In  this  section  we  want  to  consider  dispatching  segments  of  a  unit  as 
opposed  to  dispatching  a  unit's  full  capacity.  We  will  not  consider 
partial  outages  of  generating  units.  If  a  unit  is  unavailable  for  service 
then  a11_  its  segments  from  no  load  to  full  load  are  unavailable.  Likewise, 
if  a  unit  is  available  then  all  its  segments  are  available. 

In  Section  4.2,  we  found  that  in  determining  the  effective  load 
duration  curve,  it  was  only  necessary  to  consider  those  units  already 
dispatched  (to  the  left  of  the  unit  about  to  be  dispatched  under  the  load 
duration  curve).  Those  units  to  be  dispatched  afterwards  did  not  affect 
the  position  under  the  load  duration  curve  of  the  unit  about  to  be 
dispatched. 

The  same  interpretation  can  be  applied  here  with  a  slight  modification. 
That  is,  in  forming  the  effective  load  duration  curve,  any  unit  that  has 
had  at  least  one  segment  already  dispatched  can  affect  the  position  of 
the  segment  about  to  be  dispatched  and,  therefore  must  be  included.  Any 
unit  that  has  not  had  any  segments  dispatched  will  be  to  the  right  of  the 
segment  about  to  be  dispatched  and  cannot  affect  its  position.  Thus,  it 
need  not  be  included  in  forming  the  effective  load  duration  curve. 


n  "III  -  mi^^  I 


no 


Suppose  vie   are  about  to  dispatch  the  second  segment  of  unit  i.  Any 
units  that  have  had  at  least  one  segment  already  dispatched  must  be 
included  in  an  equation  of  the  form  of  equation  (4.2-14).  Unit  i  has 
already  had  one  segment  dispatched.  Should  it  be  included  in  the  effective 
load  duration  curve?  The  answer  lies  in  the  answer  to  the  following 
question.  Does  the  availability  or  unavailability  of  the  first  segment 
of  unit  i  affect  the  position  under  the  load  duration  curve  of  the  second 
segment?  While  on  first  thought  it  seems  it  does,  the  answer  is  no.  This 
is  because  if  the  first  segment  is  unavailable,  then  so  is  the  second 
segment.  Remember  that  partial  outages  are  not  being  considered.  Thus, 
if  we  are  about  to  dispatch  the  second  segment,  both  it  and  the  first 
segment  must  be  available,  and  as  a  result,  unit  i  must  not  be  included 
in  the  formation  of  the  effective  load  duration  curve. 

Most  of  the  useful  production  costing  models  to  date  (see  Booth  (1972), 
Jenkins  and  Joy  (1974),  and  Sullivan  (1974))  form  the  effective  load 
duration  curve  one  step  at  a  time  as  each  new  segment  is  dispatched.  As 
each  segment  is  about  to  be  dispatched,  that  unit  must  be  "deconvolved  out" 
of  the  effective  load  duration  curve.  This  requires  an  iterative  procedure 
that  takes  many  computations  and  may  involve  accumulating  errors  of 
the  round-off  type.  We  are  able  to  avoid  this  problem  by  forming,  at 
each  step,  the  appropriate  effective  load  duration  curve  considering  only 
those  units  that  have  had  at  least  one  segment  dispatched,  with  the  single 
exception  of  the  unit  about  to  be  dispatched.  Thus,  for  the  i-th  unit 
about  to  be  dispatched  over  its  segment  a  to  b:  (0  -  a  <  b  -  P.  „,„)» 

1  ,maX 

we  have  for  its  expected  cost 

P-^+b 

C  =  p,f^.  f   IHR.(P, -F^)ER  (P,  ),    dP,  (4.3-1) 


»eg»C»lmVStaP^  -— — .^ — ■— -; -^■••«a-^->-.J-.!?.'.— ?.-ra!su-«:.i-ii.i[,---— .^■t.r^,  ia|>f.aj/JBKi2iKtru^<M^c-^|aMHBi ii_aci;j  -  — .-miKTMunil ifc;>^JTMt  Jii •'  T-'    tm''  '— M^nrf  fi-M-^iiltrr '  U^n^ii-Fift^V  ■£!■«■  iinlln>ikc-%lir--- 


Ill 


L^  L' 


+  (q^P2P3  ...  Pi_-,)FL(PL-Pi) 

MP^q2P3  ...  Pi.,)FL(PL-P2) 
+  ... 

+  (q^qgqa  ...  qi-l)f'L^\"P^•••-^•-l)         (4.3-2) 

The  term  Pj(j=l  ,2, . . . ,i-l )  in  equation  (4.3-2)  is  the  amount  of  power  of 

the  j-th  unit  already  dispatched.  Mote  that  equation  (4.3-2)  is  M^ry 

similar  to  equation  (4.2-14),  except  P.  replaces  P.  ,  .  Also  Pt  is  now 

J         J  jmax       1 

given  by 

i-1 

Pj  =  .1  Pj  •  (4.3-3) 

J=l  ^ 

Equations  (4.3-1),  (4.3-2),  and  (4.3-3)  form  the  basic  equations  used 
in  our  production  costing  model  to  determine  expected  values  of  fuel  costs 
and  emission  rates  of  various  generating  units  when  dispatched  over  a 
particular  segment  of  their  capacity,  In  Chapter  5  we  will  apply  this 
production  costing  model  along  with  the  dispersion  model  developed  in 
Chapter  3  to  study  the  effects  on  ground  level  concentrations  of  sulfur 
dioxide  our  dispatching  scheme  can  have. 

'  4.4  Estimating  Short  Term  Maximum  Concentrations 

In  Chapter  3  we  developed  the  mathematical  model  for  estimating 
long  term  ambient  concentrations  of  sulfur  dioxide  at  ground  level  due 
to  emissions  from  electric  power  plants.  Earlier  in  this  chapter,  we 


where  EFj  (P,  )    is  now  given  by  '  j 

i-1  ■     I 
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presented  the  method  of  production  costing  and  extended  the  theory  to 
include  certain  generalized  functions  of  the  output  of  electric  generators. 
In  particular,  we  were  interested  in  long  term  estimates  of  fuel  costs 
and  emissions  of  sulfur  dioxide  as  a  function  of  the  position  of  a  segment 
of  a  generating  unit  under  the  load  duration  curve.  Since  the  load 
duration  curve  was  representative  of  some  time  period,  T,  that  is 
generally  a  month  or  longer,  these  estimates  or  expected  values  were 
considered  to  be  long  term  estimates. 

Before  discussing  the  application  of  the  models  of  Chapters  3  and 
4  to  a  realistic  size  electric  utility  system  (Chapter  5),  we  want  to 
present  and  discuss,  in  this  section,  a  method  for  predicting  maximum 
short  term  concentrations  (on  the  order  of  hours)  from  the  estimates  of 
long  term  average  concentrations  of  sulfur  dioxide.  The  method  was 
developed  by  Larsen  (1969).  We  will  make  use  of  his  method  to  extend  the 
usefulness  of  our  model  in  predicting  both  long  term  and  short  term 
estimates  of  ground  level  concentrations  of  sulfur  dioxide. 

In  Chapter  1  we  noted  that  the  U.S.  Government  through  the  EPA  has 
issued  both  long  term  standards  (annual  averages)  and  short  term  standards 
(3  hour  and  24  hour  maximum  concentrations).  Although  we  will  be  primarily 
concerned  with  the  long  term  averages,  this  method  allows  us  also  to  report 
the  predicted  short  term  maximum  concentrations.  Thus,  we  could  elect  to 
control  on  either  the  short  term  or  long  term  concentrations  (Chapter  5). 

In  this  section  we  will  develop  those  aspects  of  Larsen's  method  that 
we  find  of  use  in  our  model.  The  complete  development  and  detailed 
derivation  is  given  by  Larsen  (1971).  Since  Larsen's  method  assumes  the 
lognormal  distribution,  the  reader  is  referred  to  Appendix  A  for  a  dis- 
cussion of  the  mathematics  of  the  lognormal  distribution  as  well  as  seme 
of  its  properties. 


iini^'iiiB'H'M —    I '  li'ii  riM  ifn     ili^iB  I'l    —     .  ^  •.--.j,^      i^iiB    ■■■■-||-     r-trm'a    tn^MBJiltun'i  if.lPii'iri»^MiiinHiii^^~-iinrrhli-(liiriliii7TilH-T~* 
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Larsen  (1971)  in  reference  to  his  method  states,  "The  two  dominant 
features  of  this  new  method  are  derived  from  observations  that  indicate 
(1)  air  pollutant  concentrations  are  approximately  lognormally  distributed 
for  all  pollutants  in  all  cities  for  all  averaging  times;  and  (2)  median 
concentration  is  proportional  to  averaging  time  to  an  exponent."  (p.  1) 
The  observations  referred  to  above  were  conducted  by  recording  continuous 
air  pollutant  concentration  data  for  seven  pollutants  (including  sulfur 
dioxide)  in  six  cities  over  a  three-year  period  (see  Larsen  et  a1 .  (1967)). 

During  this  three-year  period,  Larsen  et  al .  (1967)  performed  a  study 
which  showed  that  concentrations  of  various  pollutants  follow  a  lognormal 
distribution.  If  one  plots  the  logarithm  of  the  concentrations.  In  x> 
versus  the  number  of  standard  deviations,  z,  from  the  median,  then  a  straight 
line  results  in  the  form 

In  X  -  oz+u 
where  a   is  the  standard  deviation  and   is  the  mean  of  the  normal 
distribution  after  the  transformation  (A-1)  from  Appendix  A.  Thus,  a  is 
seen  to  be  the  slope  of  the  line,  and  y  is  the  intercept.  Using  equations 
(A-7)  and  (A-11),  this  becomes 


lnx=2lns+lnm 

3  y 


or 


X        =      "IgSg^  (4.4-1) 

where  m  is  the  geometric  mean,  and  s  is  the  geometric  standard  deviation 
(see  Appendix  A). 

Equation  (4.4-1)  is  valid  only  when  all  quantities  are  related  to  the 
same  averaging  time.  Thus  for  24  hour  concentrations,  we  must  have  24  hour 
values  for  m  and  s  .  '^'hen  we  have 

g    g 

^■24  ^  ''g24  '^'g24'^  ^4.4-2) 


':>«c»niWiat4i««i;»u^an^9'Savis7=''rtf''«ni'?dfl7llli£«i;aUEti  t  irrniim  liTniiii  <-ir>Mwirr^^°^'fii^-Mr»i  ■--^-'^-^^^■--'-  -  ~^:  .-...■- ^--    .    ^  .i^sniikJt  i        n       il   i     1 1  ■^■im    r       -  i  ■  ■  ■  ■■!■■  ■  ii  ii  -ibih    h  i  -     ■■■hhmmt-thii  "nn-  ar^- 
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Since  we  will  be  interested  in  various  short  term  averaging  times, 
it  will  be  necessary  to  convert  m  and  s  from  one  averaging  time  to 
another.  Two  characteristics  on  which  the  method  is  based  are  stated  here: 

1)  Median  concentrations  are  proportional  to  averaging  time  to 
an  exponent. 

2)  For  the  longest  time  period,  t   ,  the  arithmetic  mean  equals 
the  median. 

Since,  as  is  shown  in  Appendix  A  for  lognormal  distributions,  the 
median  equals  the  geometric  mean,  we  can  express  characteristic  one,  above,  as 
m.  =  kt^  (4.4-3) 

Taking  logarithms  of  both  sides  of  equation  (4.4-3)  we  have 

In  m  =  In  k  +  P  In  t  (4.4-4) 

g 

Equation  (4.4-4)  would  plot  as  a  straight  line  on  log-log  paper.  Choosing 

three  points  on  that  line  representing  the  total  averaging  time,  t  ., 

averaging  time  t  ,  and  t,  ,  we  can  write  three  equations.  These  are 

In  m  tot  =  In  k  +  P  In  t^  ^  (4.4-5) 

g  tot 

In  m   =  In  k  +  P  In  t  (4.4-5a) 

ga  a 

In  m  .  =  In  k  +  P  In  t^^  (4.4-5b) 

Using  characteristic  two  the  first  of  these  three  may  be  written  as 

In  m  =  In  k  +  P  In  t^^^  (4.4-5) 

where  m  is  the  mean  concentration  for  the  total  time  period. 
Using  equation  (4.4-5)  and  (4.4-5a)  we  can  write 

ln(m/mg^)  =  P  In  (t^^^/t^) 

or 


cSBTtr  -m'  BWg=i-aHn'aig<r'Sf=:^a'!i^aSfiBr.'«~»»'V^>k*'gRja  -. .. — - . . ..— -.. — .<=j--j,j:.^  ^^.~,r.^  ^^^,^^^j^,.  .^^  -.-  -r=-  ^  -—.-  l-  --    -h  ^i^  twrr.  Mflifi-riii  i  v^,  m  -'-i»ii  --T  r  r  'T'h'  7  •  r    — i  i  »ir  n  t  —  -■    r  1 1  irr '  v  nn  iff  -  r  "*n  ntT~r~'  r'r 


IIS 


ln(rn/m      ) 
'"''Hot'   a^ 

Similarly,   choosing  equation   (4.4-5b)  we  can  write 

ln(m/m  ,  ) 
P  =  ln(t       /t  )  ("^-^8) 

'"^Hot^  b^ 

Using  equation  (A-12)  »  which  is  repeated  here  for  convenience, 


I 


we  have  i 

12-  [ 

E[x]  =  exp(y+  |a^)  (At12) 

Then  since  the  mean  of  the  concentrations  is  also  the  expected  value, 
we  have 

m  =  E[x] 

also  recall  that 
u  =  In  m 

g 

and    a  =  In  s 

g 

Taking  logarithms  and  using  the  above  substitutions,  we  can  write 

1       2 
In  m  =  In  m  +  2  (In  s  )  (4.4-9) 

Using  equation  (4.4-9)  and  equating  equations  (4.4-7)  and  (4.4-8) 
yields 

(in  Sg,)'  "  ^"'Wb'  I 

I' 
f 

or  finally  ! 

V  =  (^gb^"  (4.4-10) 

i 

where 


gywffltWIPfrSII-^^Ttf^tt*?^^.^  ■^M  ili  ^B-JUlff?a''li'IW-ftri«^iet»'ZTW£Kf>«^g^Ma<i<'»-»^*^  1liMji.^j*gyTi>.,^y*«ifjwy>. 
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Equations  (4.4-10)  and  (4. 4-11)  relate  the  geometric  standard  deviation, 

s_.,  for  averaging  time,  t  ,  to  the  geom.etric  standard  deviation,  s  ,  , 
9°  a  gb 

for  averaging  time,  t^^.  Also  t    is  the  total  averaging  tine  period 

(usually  a  year).  For  clarity,  we  will  present  an  example  involving  all 

these  quantities  at  the  end  of  this  section. 

One  can  also  develop  a  relationship  between  the  geometric  means, 

similar  to  equations  (4.4-10)  and  (4.4-11).  However,  it  is  not  necessary 

to  do  so,  since  equation  (4.4-9)  relates  the  geometric  mean  to  the  geometric 

standard  deviation.  Thus,  we  need  only  change  the  time  base  for  s  and  ■ 

then  relate  m  to  it,  through  equation  (4.4-9).  That  is,  rewriting 
9 

equation  (4.4-9),  we  have 

%  =  —7177 — ^zr  (^•4-12) 

^^       exp[2-(ln  Sgg)  ] 

The  use  of  equation  (4.4-1)  for  predicting  short  term  concentrations 
depends  on  z,  the  number  of  standard  deviations  from  the  median  that  the 
concentration  is  located.  This  distance,  in  turn,  depends  on  the  relative 
frequency,  f^,  with  which  the  particular  concentration  is  expected  to 
occur.  For  example,  the  highest  one  hour  concentration  in  one  year 
(8,760  hours)  would  occur  with  a  frequency  of 

For  ranking  the  order  of  highest  concentrations,  Pearson  and 
Hartley  (1962)  have  suggested  the  following  formula 

^,  -  ^  (4.4-1.3) 


itivi«iNli!ii>;^^n«na»'afJg»Eiii;-ii>i  *<.»^^v»J:^^gu^^<tg^Btfn^a*Jl■rt■^rl^'^^^■^l^liIPMJifg^lll!l■ac^g>.=^ 


117 


where  f  is  the  frequency  (as  a  decimal),  n  is  the  number  of  samples, 
and  r  is  the  rank,  e.g.,  r=l  is  the  highest,  r=2  is  the  second  highest, 
etc.  For  our  purposes  equation  (5.2-13)  becomes 

t 

f  =  -^-   (r-.4)  (4.4-14) 

^   ^tot 

where  t  and  t  ,  (in  hours)  are  as  described  earlier, 
a     tot 

The  number  of  standard  deviations  from  the  mean  can  then  be  written 


as 


l-f^ 


z       1  2 
exp  (-  2"  2  )  "^s 


or 


0.5-f  =  I  exp  (-  1  s^)  ds  (4.4-15) 

•'o 

We  have  solved  equation  (4.4-15)  for  z,  for  various  values  of  t 

(which  relate  to  f  by  equation  (4.4-14))  with  t.^  equal  to  8,760  hours 

(1  year).  Table  3  shows  the  z  values  for  the  first  and  second  highest 

concentrations  for  several  values  of  t  (in  hours).  The  solution  of 

a 

equation  (4.4-15)  is  incorporated  into  the  computer  program  discussed  in 
Chapter  5,  so  that  highest  concentrations  for  any  averaging  time  of 
interest  may  be  predicted. 

\'ie  will  finish  this  chapter  with  two  examples  that  illustrate  the 
use  of  the  pertinent  equations  of  this  chapter,  as  they  are  used  in  the 
computer  algorithm  of  Chapter  5. 

Example.  (1 ) 

Suppose  we  are  given 

,    3      ,  .« 

m  =  45  (ug/m  annual ) 

and 


-iTi~t~"~y'i — rir"iiiii'ri'tiiiTTr  i  m  rr  "i  if-r  ~>i'i"t  ■  iiti  n'l  ■  'H'^iriiiifr  -iiTrini — i   in  '~Ti'^   i»r  i-ii  mi     ir    — *    i  T""!'  rn  ir  n  mTtr-  ^'rc    •to-»-=n--"-*-*^i-**'^«Mi>'i'^i,iiBi«rt>»>*rir*i  "-Tir*^"*^!^ 
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Table  3,     Number  of  Standard  Deviations  from  the 
Median. 


t     in  hours  z1 ,   1-st 

highest 


1  3.81 

2  3.64 

3  3.53 

4  3.46 


z2,   2-nd 

highest 

3.57 

3.42 

3.26 

3.18 

''V.mem'tt  ii iwi'"  ^"'  ■  > * ' f  h ^t*! ' pfc*>i^t-i»iiraiii  jf^ttrj^msr^-n  ■-■  Miiaw>ffiftifiir^s>i<!gi^irt':MiiL wp«e~iawii *Migiwi ■■  itf"*»' ti*  n * •i«ifHgrit..»JK3P'«3m "iiwy  iaii»rt >w/M.i*qw*af«f»iitf*rriMi»«»i^wii4ifi 
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Sg24  =  2.21 

and  we  want  to  determine  the  highest  and  second  highest  three-hour 
concentrations.  Since  the  arithmetic  mean,  m,  is  an  annual  value,  the 
total  time  period,  t  . ,  is  8,760  hours.  Note  that  the  arithmetic  mean 
is  the  output  of  our  long  term  dispersion  model  (Chapter  3), 

The  first  step  is  to  find  s^  by  use  of  equations  (4,4-10)  and 
(4.4-11).  This  is  done  as  follows 

s   =  (s   )" 

where 

n  -  rln  (8,760/3)  ,1/2  _  ,  . 
"  -  Mn  (8,760/24)J    -  '-^^ 

Then 

s  ,  =  (2.21)^-^^  =  2.51 

g3 

m 

and  from  equation  (4.4-12) 

m  -  =  T — ^ ^  =  29.41 

9^   exp[^(ln  2.51)^] 

To  determine  z  we  use  Table  3  and  find 

z  (highest)  =  zl  =  3.53 
and 

z  (second  highest)  =  z2  =  3.26 

Then  from  equation  (4.4-1) 

X3  (highest)  =  mg3  {s^^f^ 

X3  (highest)  -  29.41  (2.51)^-^^ 
or 


MMilMitefii^Jwpiia  m,  I*    igB7|iKPir^*^wiiii[!PS1f'w>:;!T=''J»tJirii'»^-rslw»l»i.-^«itfT^ 
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X3  (highest)  -  752.7  (yg/m^) 

3 
.which  compares  favorably  with  the  1,300  (yg/m  )  three-hour  maximum 

standard  not  to  be  exceeded  more  than  once  per  year. 

For  the  second  highest  we  have 


X3  (second  highest)  =  29.41  (2.51)^-^^ 

q 

X-j  (second  highest)  =  594.6  (^g/m  ) 


Example.  (2) 

As  mentioned  above.  Table  3  is  based  on  a  total  averaging  time  of 
one  year.  To  illustrate  its  use  for  other  total  averaging  times,  suppose 
the  arithmetic  mean  in  Example  (1)  is  for  a  three-month  period  (2,190 
hours).  Further,  suppose  tliat  we  desire  estimates  of  the  highest  and 
second  highest  one-hour  concentrations.  Then  we  have 

_  pin  (2,190/1)  1/2 
"  "  Hn  (2,190/24)^    ""  '•■^' 


Then 


Sg,={Sg2,)"=  (2.21)^-2' 


S  T  ~  2 . 82 


and 


45 
m  ,  =  ^ -K-    =  26.33 

^  .   exp  ij   (In  2.82)^] 

Notice  that  the  ratio  of  one  hour  to  2,190  hours  is  the  same  as 
four  hours  to  8,760  hours  (one  year).  Than  using  Table  3  for  four  hours, 
we  have 

z  (highest)  =  zl  =  3.46 

z  (second  highest)  =  z2  =  3.18 


-r»»iu»in.tMiTTtt*-rt  ■ 
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rhus 


X3  (highest)  =  26.33  (2.82) 
<3 


7   0,^3.46 


Xo  (highest)  =  945.3  (  g/m^) 


while 

X3  (second  highest)  =  26.33  (2.82)-^-^^ 

X3  (second  highest)  =  708.2  (  g/rr?) 

This  completes  this  chapter.  In  the  next  chapter  we  put  the  models 
of  Chapters  3  and  4  together  into  a  computer  program  to  simulate  a 
realistic  size  electric  power  system. 


CHAPTER  5 
DISCUSSION  OF  THE  DIGITAL  COMPUTER  PROGRAM 

In  this  the  final  chapter,  we  will  discuss  the  computer  program 
we  have  developed  for  selecting  the  optimal  strategy  for  committing 
electric  generating  units  so  as  to  meet  the  Federal  Ambient  Air  Quality 
Standards.  As  pointed  out  in  the  earlier  chapters,  this  approach  involves: 
1)  estimating  the  ground  level  concentrations  of  sulfur  dioxide  from 
electric  power  plants  through  a  steady  state  long  term  dispersion  model 
(Chapter  3),  2)  selecting  an  optimal  unit  commitment  strategy  through 
the  extended  method  of  production  costing  (Chapter  4),  and  3)  estimating 
the  highest  and  second  highest  concentrations  via  Larsen's  technique 
(Chapter  4).  The  actual  technique  for  selecting  the  optimal  commitment 
strategy  has  been  purposely  delayed  until  this  chapter.  In  the  next 
section  this  technique  is  discussed,  along  with  the  other  salient  features 
of  the  computer  program. 

5.1  The  Digital  Computer  Program 

The  most  expedient  method  for  discussing  the  computer  program  used 
in  this  work  is  to  follow  the  sequence  of  steps  executed  by  the  computer 
and  outlined  in  Figure  15.  The  flow  diagram  (Figure  15)  begins  with  the 
input  data  and  ends  with  the  output  results.  Both  of  these  are  discussed 
thoroughly  in  Section  5.2.  For  these  reasons  the  discussion  here  begins         j 
with  the  second  block  in  the  flow  diagram--the  computation  of  the 
environmental  coefficients. 
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Figure  1  5.  Flow  Diagram  of  Computer  Program 
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The  environmental  coefficients  are  also  called  influence  coefficients 
and  sometimes  called  "x  over  Q"  coefficients  (Turner  (1970)).  These 
definitions  come  from  equations  (3.4-12)  and  (3.4-13).  If  x  is  divided  by 
Q  the  resulting  equation  on  the  right  is  only  a  function  of  the  meteorological 
parameters  and  the  effective  stack  height.  This  is  shown  below  for  equation 
(3.4-13) 

X  =  2^  exp  [-  1  (f  )2]  (5.1-1) 

Q   a. Xu        ^     ^z 

The  term  on  the  right  of  equation  (5.1-1)  is  independent  of  the 
emissions,  Q.  The  concentration,  due  to  a  particular  generating  unit, 
can  be  determined  whenever  both  the  environmental  coefficient  and  the 
emissions  of  the  particular  unit  are  known,  thus  separating  the  dispersion 
calculations  from  the  production  costing  calculations.  The  dispersion 
calculations  need  only  be  made  once  for  a  given  set  of  meteorological 
conditions.  As  soon  as  the  input  data  is  complete,  the  program  computes 
an  environmental  coefficient  for  each  grid  (see  Section  5.2)  due  to  each 
generating  unit  and  stores  these  results  for  future  use. 

Next,  the  control  is  set  to  economics,  and  the  program  selects  the 
unit  commitments.  This  is  done  by  committing  the  unit  segments  to  serve 
the  load  (Chapter  4)  in  an  order  that  minimizes  the  total  system  fuel 
costs.  That  is,  at  each  step  each  unit  in  a  group  (Section  5.2)  is  a 
candidate  for  supplying  the  next  increment  (segment)  of  energy.  Equation 
(  4.3-1)  is  used  in  the  program  to  compute  the  fuel  costs  of  each 
candidate  (unit),  and  the  one  whose  fuel  cost  is  the  minimum  is  selected 
and  placed  in  that  position  under  the  effective  load  duration  curve. 
This  process  is  continued  until  all  generating  units  have  been  committed. 


Obviously,  as  the  last  segment  of  a  unit  is  committed,  it  is  no  longer  a 
candidate  for  any  further  positions.  Since  the  criterion  for  selecting 
a  unit's  segment  for  position  under  the  effective  load  duration  curve 
is  based  on  fuel  costs,  the  control  is  said  to  be  based  on  economics. 

The  order  in  which  the  units  are  conmitted,  as  well  as  the  individual 
segment  fuel  costs,  are  stored  for  future  reference.  The  individual 
segment  emissions  can  then  be  determined  from  equation  (4.1-21),  which 
is  repeated  here  for  convenience 

Q  =  (KfJC  ■  (4.1-21) 

These  emissions  are  used,  along  with  the  environmental  coefficients,  in 
the  next  block  of  the  flow  diagram  to  check  for  air  quality  violations. 
A  violation  is  considered  to  have  occurred  if  any  grid  has  a  concentration 
level  is  excess  of  60  percent  of  the  annual  primary  standard  of  80 
micrograms  per  cubic  meter.  The  selection  of  60  percent  is  in  keeping 
with  a  conservative  strategy;  and  the  actual  percentage  level  could  easily 
be  changed  to  any  other  level  desired.  If  no  violations  have  occurred, 
the  program  outputs  the  results  and  terminates. 

If  at  least  one  violation  occurs,  the  control  is  set  to  air  quality 
as  is  shown  in  the  flow  diagram;  and  the  program  again  selects  unit 
commitments.  However,  in  this  case,  the  selection  is  based  on  air  quality 
in  the  following  manner.  For  each  candidate  for  a  given  position  under  the 
effective  load  duration  curve,  the  program  computes  the  resulting  emissions 
of  sulfur  dioxide.  Next  using  the  environmental  coefficients,  the  program 
computes  and  sums  the  ground  level  concentrations  of  sulfur  dioxide  that 
would  occur  throughout  all  the  grids  (geographical  area)  for  each  candidate 
occupying  that  position  under  the  effective  load  duration  curve.  That 
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unit's  segment  which  produces  the  minimum  total  ground  level  concen- 
trations is  selected.  This  process  is  continued  until  all  the  units  have 
been  committed,  thus  defining  the  desired  commitment  order  and  associated 
segment  fuel  costs  for  future  use.  Since  the  individual  emissions  can 
easily  be  recovered  through  the  use  of  equation  (4.1-21),  it  is 
unnecessary  to  store  them.  The  selection  process  is  based  on  minimizing 
total  concentrations  of  sulfur  dioxide;  hence,  the  control  is  said  to  be 
based  on  air  quality. 

Next,  the  program  checks  again  for  air  quality  violations.  If  there 
are  none,  the  control  of  the  program  is  routed  to  the  block,  "search  for 
best  case."  If  there  are  violations,  the  control  is  routed  to  the  block, 
"compute  new  sulfur  contents."  We  will  discuss  both  these  blocks  beginning 
with  the  latter  block  first. 

If  there  is  at  least  one  violation,  there  is  no  way  to  avoid  air  quality 
violations  by  reordering  the  unit  commitments,  since  all  units  have  been 
committed  in  an  order  that  minimizes  the  total  concentrations.  Thus,  if 
the  energy  is  to  be  produced,  the  emissions,  and  hence  the  sulfur  contents 
must  be  lowered.  There  are  various  schemes  for  deciding  which  units  should 
have  their  sulfur  contents  lowered;  for  example,  a  Lagrange  multiplier 
formulation,  that  minimizes  increased  fuel  costs  due  to  more  expensive  low 
sulfur  fuels,  subject  to  the  inequality  that  the  grid  concentrations  all 
be  less  than  or  equal  to  the  60  percent  standard.  While  this  method  is 
preferred  in  theory,  it  is  of  little  practical  use  here.  This  is  because 
in  practice  only  those  units  that  do  not  suffer  significant  deratings 
(Chapter  2)  are  actually  candidates  for  lower  sulfur  fuels.  The  decision 
as  to  which  units  must  burn  lower  sulfur  fuels  is  likely  to  be  made  on  the 
basis  of  avoiding  combustion  problems  rather  than  on   increased  fuel  costs. 
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For  this  reason,  we  have  selected  the  procedure,  favored  by  the  EPA 
(1975),  of  a  direct  roll-back  of  emissions  on  all  sources  capable  of  burning 
lower  sulfur  fuels.  As  discussed  in  the  next  section,  the  user  can  specify 
whether  a  generating  unit  is  capable  of  lower  sulfur  contents  or  not,  as 
well  as  a  minimum  low  sulfur  content.  The  program  assumes  a  three  percent 
increase  in  fuel  prices  for  each  one  quarter  percent  reduction  in  sulfur 
content.  This  amounts  to  a  16  percent  increase  in  fuel  costs  for  a  full 
percent  reduction  in  sulfur  content  (PEDCo  (1975)). 

In  the  program  the  procedure  is  to  select  the  grid  having  the  highest 
concentration  (most  serious  violation)  and  to  roll -back  the  emissions  of 
the  preselected  units  to  bring  this  grid  into  conformance  with  the  60  percent 
standard.  Next,  the  new  sulfur  contents  are  lowered  further  to  the  nearest 
one  quarter  percent.  For  example,  3.26  percent  is  lowered  to  3.25  percent; 
while  3.24  percent  is  lowered  to  3.00  percent  sulfur  content.  This  is  in 
keeping  with  a  conservative  approach  and  also  to  avoid  specifying 
unrealistic  numbers,  such  as  3.17  percent. 

Since  it  is  possible  that  some  grid  violations  are  caused  by  units  not 
selected  for  lower  sulfur  contents,  the  program  checks  to  insure  that  ari_ 
violations  have  been  brought  into  conformance.  If  some  violations  are 
still  found  to  be  present,  the  procedure  outlined  in  the  paragraph  above  is 
repeated  until  all  violations  have  disappeared.  Also,  when  the  program 
control  is  first  transferred  to  "compute  new  sulfur  contents,"  a  check  is 
made  to  see  if  the  violations  can  be  removed  by  placing  all  selected  units 
at  their  specified  minimum  low  sulfur  content.   (If  not  specified,  a  0.5 
percent  minimum  is  assumed.)  If  violations  still  occur  at  this  point,  the 
program  terminates  and  informs  the  user  that  it  is  not  possible  to  meet 
the  60  percent  standard  under  the  conditions  he  has  specified. 
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Notice  in  the  flow  diagram  that  after  computing  the  low  sulfur  con-         i 
tents,  the  control  is  transferred  to  the  block  "set  control  to  economics."       | 

r 

This  is  done  since  we  are  assured  that  no  violations  will  occur  with  all         : 
units  committed  on  the  basis  of  air  quality,  or  else  the  program  would  have      i 


terminated.  Since  the  lower  sulfur  contents  have  been  lowered  to  the 
nearest  quarter  of  a  percent,  it  is  possible  that  no  violations  will  occur 
with  units  committed  on  the  basis  of  economics.  This  is  checked;  and  if 
none  are  found  to  occur,  the  program  outputs  the  results,  which  includes 
the  commitment  order  based  on  economics  and  the  new  lower  sulfur  contents. 

If  at  least  one  violation  is  found  to  occur,  the  program  transfers 
control  to  the  block,  "search  for  best  case,"  which  is  discussed  next. 
Notice  that  the  flow  diagram  is  drawn  so  that  the  program  reselects  the 
order  based  on  air  quality.  This  is  because  changing  the  sulfur  contents 
may  change  the  commitment  order.  Actually,  this  repeated  selection  is  per- 
formed while  in  the  block  "compute  new  sulfur  contents,"  to  insure  that 
there  are  no  violations  occuring  based  on  both  the  new  lower  sulfur 
contents  and  the  new  commitment  order. 

The  final  block  to  be  discussed  in  the  flow  diagram  in  this  section 
is  the  block,  "search  for  best  case."  This  block  is  necessary  since  when 
control  is  transferred  to  this  block,  two  important  facts  are  known: 

1)  Violations  will  occur  if  units  are  committed  on  the  basis  of 
economics. 

2)  Violations  will  not  occur  if  units  are  committed  on  the 
basis  of  air  quality. 

As  discussed  in  the  next  section,  the  units  are  placed  in  groups,  and 
each  group  of  units  is  either  committed  on  the  basis  of  economics  or 
air  quality.  If  there  are  three  •jroups,  there  are  eight  (2")  cases 
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(possibilities)  to  be  considered.  The  first  case  (all  on  economics)  ; 

i 

has  been  considered,  and  the  eighth  case  (all  on  air  quality)  has  also  t 

1 
been  considered.  Since  the  program  has  stored  the  conmitment  orders  of         [ 

both  these  cases,  it  is  possible  to  consider  the  other  six  cases  by  [ 

I 

"swapping"  in  and  out  the  commitment  order  of  the  units  within  groups.  I 

i 

This  is  done  in  the  program.  ; 

All  eight  possible  combinations  of  the  three  groups  being  committed         , 

I 

on  economics  or  air  quality  are  considered  and  air  quality  violations  are 
checked.  For  those  cases  that  do  not  have  any  air  quality  violations,  the 
fuel  costs  are  computed  and  the  "best  case"  is  selected  as  the  one  with  the 
minimum  fuel  cost.  If  all  the  cases  are  found  to  contain  air  quality 
violations,  then  the  case  selected  is  the  one  in  which  all  groups  are 
committed  based  on  air  quality  (case  8  for  this  example),  and  control 
passes  to  the  block,  "output  the  results." 

Before  discussing  the  input  and  output  blocks  in  the  next  section, 
we  need  to  consider  the  dispersion  model  near  the  generating  source.  As 
mentioned  in  Section  3.4,  it  must  be  modified  for  distances  closer  than  one 
kilometer  to  the  source.  We  find  it  more  convenient  to  make  this 
modification  here,  since  it  goes  well  with  the  discussion  of  the  computer 
program. 

In  the  program,  sources  are  placed  at  the  center  of  a  square  grid  of 
dimensions  It   by  It.     For  example,  if  r  equals  2500  (m) ,  then  the  grid 
is  5000  (m)  on  each  side.  We  will  consider  the  grid  to  be  roughly  a  circle 
of  radius  r,  for  the  purposes  of  obtaining  the  average  concentration 
throughout  the  grid  due  to  the  source  in  the  grid.  This  source  square 

average  concentration  will  be  given  by 

16    1  ri^x(x,f.) 

Xsc  ^  )  ^  ~   !  f~-dx  (5.1-2) 

^-   i  =  l   '  "^  JO    1- 
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where  x(x-,f^)  is  given  by  equation  (3,4-12)  or  (3.4-13),  r  is  one  half 
the  grid  linear  dimension,  and  f.  is  the  frequency  with  which  the  wind 
blov/s  tov/ard  the  i-th  sector.  Note  that  the  f.  in  equation  (3.4-12)  or 
(3.4-13)  identifies  the  particular  term  as  representing  the  i-th  sector. 
Equation  (5.1-12)  can  be  rewritten  as 

X_  =M  X(x)dx  il     f  ) 
ss   r  jQ        .^^  1 

or 

1  T- 

since  the  summation  equals  one,  and 

X(x,f.) 
X(x)  =  ~f7-^  (5.1-4) 

For  the  case  of  no  inversion  layer  aloft  equation  (5.1-4)  becomes 

X(x)  =%^  exp  [-^(-^)2]  (5.1-5) 

Notice  that  in  equation  (5.1-2)  it  was  only  necessary  to  integrate 
over  the  x-coordinate.  The  integration  over  the  y-coordinate  was  done 
in  Chapter  3.  Also  each  f.  is  treated  as  the  probability  of  concentrations 
appearing  in  the  i-th  sector  due  to  the  source  in  the  source  square. 

Hanna  (1972)  derives  the  source  square  long  term  concentration 
formula  as 


^  .i,2i,  r^ 


ss  _.Z   (f?)  J„  '<(^'<''<  (5-1-6) 


Ttr      •'0 

where  the  details  can  be  found  in  Henna  (1972).  The  ur     tenri  is  due  to 
the  area  of  the  area  of  the  circle  of  radius  r.  Notice  that  the  units  of 
equation  (3.1-6)  are  yg/m  ,  whicli  appears ' incorrect.  Perhaps  the  discrepancy 


£.^£;b?e-.-^-&iii.~i. 
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'An  unpublished  report  available  from  the  author. 


is  due  to  summing  over  the  area  and  then  dividing  by  the  area,  when  the 
crosswind  distance  has  been  summed  and  divided  already  in  the  averaging 
process  of  Chapter  3.  We  believe  it  is  only  necessary  to  sum  and  divide        j 

by  the  linear  distance  r  as  we  have  done  in  equation  (5.1-3).  Notice  [ 

3  ■ 

also,  the  units  in  equation  (5.1-3)  are   correct  (yg/m  ).  I 

5 . 2  Discussion  of  the  Application  of  the  Computer  Program  to 

a  Realistic  Size  System  I 

i 

The  digital  computer  program  discussed  in  this  chapter  has  been  \ 

written  as  a  source  program  in  the  BASIC  language.  It  is  an  interactive 

program  in  the  sense  that  the  program  prompts  the  user  for  the  input  data  [ 

and  asks  the  user  to  make  decisions  about  whether  certain  less  relevant  j 

output  should  be  printed.  The  program  has  been  operating  on  a  Hewlett  i 

Packard  Access  System  2000  time-shared  mini-computer,  but  should  be 

compatible  with  any  digital  computer  supporting  the  BASIC  language  with 

i 
file  structures.  For  the  following  discussion,  the  reader  is  referred  to 

the  report.  User's  Guide  to  ENVECO"'",  where  a  detailed  description  of  the 

use  of  the  program  is  given. 

The  program  begins  execution  by  reading  the  input  data  from  a  file 

(BASIC  formatted  file).  The  user  is  then  given  the  opportunity  to  change 

any  input  data.  Any  changes  to  the  input  data  are  automatically  recorded 

in  the  file.  In  this  manner  each  successive  use  of  the  program  has  the 

user's  most  recent  configuration.  The  details  of  constructing  the  user's 

initial  configuration  (also  called  base  case)  are  given  in  the  User's 

Guide  to  EMVECQ.  The  program  has  been  written  with  the  flexibility  to 
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return  to  the  user's  base  case  at  any  time.  This  feature  allows  the  user 

a  simple  procedure  for  undoing  configuration  changes  in  the  input  data  | 

he  mav  no  longer  desire.  1 

i 

As  mentioned  above  as  soon  as  the  program  reads  in  the  input  data,  [ 

the  user  is  given  the  opportunity  to  make  changes  in  the  input  data.  The  [ 

procedure  for  doing  this  is  discussed  next.  The  input  data  are  divided  [ 

into  two  categories,  system  data  and  generating  unit  data.  We  will  discuss  [ 

the  system  data  first.  I 

I 

The  system  data  are  further  divided  into  six  sections.  The  first  ' 

section  is  the  load  duration  curve  data.  These  data  represent  100  ordinate      | 

values  equally  spaced  along  the  abscissa  from  zero  to  the  maximum  expected       I 

i 
system  load,  P     .  The  program  asks  the  user  if  he  wishes  to  see  or         I 

L,maX  ■ 

1 

change  any  of  these  data.  If  so,  he  can  do  so;  if  not  the  program  continues      j 
to  the  next  section,  the  wind  rose  frequencies. 

Beginning  with  this  second  section,  each  remaining  section  of  system 
data  are  displayed  at  the  user's  terminal  and  the  user  is  asked  for  any 
changes.  All  questions  asked  by  the  program  should  be  answered  yes  or  no. 
The  program  is  written  such  that  any  answer  other  than  "no"  is  automatically 
interpreted  as  "yes".  This  means  that  in  response  to  a  question  on  whether 
certain  data  should  be  printed,  a  mistaken  entry  (other  than  no)  would 
default  to  yes  and  print  out  the  data.  Also  in  the  input  stream,  a  mistaken 
enter  (other  than  no)  leaves  the  user  with  the  opportunity  to  enter  the 
new  data  rather  than  have  the  program  continue  executing  with  the  wrong 
(old)  data. 

In  the  second  section  of  system  data,  the  program  prints  the  16 
wind  rose  frequencies  from  north-northeajt  to  north.  The  user  should 
realize  that  a  change  in  cne  frequency  may  chcnge  other  frequencies.  The 
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program  automatically  sums  the  16  frequencies  and  divides  each  one  by 
the  sum  to  insure  that  the  sum  of  the  16  wind  rose  frequencies  is  equal 
to  one. 

In  the  third  section  of  system  data,  the  program  inquires  as  to 
whether  or  not  there  is  a  stable  layer  aloft.  If  not,  the  mixing  depth 
is  considered  infinite  (see  Chapter  3)  and  equation  (3.4-13)  is  selected. 
If  yes,  equation  (3.4-12)  is  selected,  and  the  program  prompts  the  user  to 
enter  the  height  of  the  stable  layer  in  meters. 

In  the  fourth  section  of  system  data,  the  program  prints  the  time 
period  in  hours  for  which  the  simulation  is  to  be  representative,  the  grid 
size  in  meters,  and  the  number  of  grid  rows  and  columns.  The  geographical 
area  is  considered  to  consist  of  rows  and  columns  of  square  grids.  The 
maximum  number  of  12  rows  and  12  columns  could  be  changed  in  the  dimension 
statement.  Each  source  (generating  unit)  is  considered  to  reside  at  the 
center  of  its  grid;  and  each  grid  concentration  is  computed  for  the  center 
of  the  grid  and  considered  representative  throughout  that  grid.  A  maximum 
of  16  generating  units  can  be  accommodated.  After  viewing  these  data, 
the  user  may  change  any  or  all  of  them  (or  none  of  them). 

In  the  fifth  section  of  the  system  data,  the  program  prints  the  10 
meter  wind  speed  in  meters  per  second,  the  ambient  temperature  in  degrees 
Kelvin,  the  number  of  groups  of  units,  and  the  number  of  segments  the  units 
have  to  be  committed  (Chapter  4),  A  maximum  of  four  groups  of  units  is 
possible  and  each  group  can  contain  up  to  four  generating  units.  The 
placing  of  units  into  a  group  to  be  com.mitted  before  any  units  are  committed 
in  the  next  group,  is  a  trade-off  between  accuracy  and  computation  time. 
Inspection  of  equation  (4.3-2)  shows  that  the  effective  load  duration  curve 
becomes  unwieldy  for  more  than  about  four  units- 
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In  the  sixth  and  final  section  of  the  system  data,  the  program  prints 
the  number  of  generating  units  in  each  group.  As  in  all  the  other  sections, 
the  user  can  elect  to  change  these  numbers.  This  completes  the  system  data. 

The  second  category  of  input  data,  the  generating  unit  data,  is 
divided  into  three  sections  for  each  generating  unit.  This  category  begins 
with  the  program  asking  the  user  if  he  wishes  to  change  any  of  the  generating 
unit  data.  If  the  answer  is  no,  then  the  computations  begin  (see  Section 
5.1).  If  the  answer  is  yes,  the  program  asks  for  the  group  number  and  unit 
number  of  the  generating  unit,  and  displays  the  current  data  for  that  unit 
in  three  sections. 

In  the  first  section  of  generating  unit  data,  the  program  prints  the 
maximum  output  in  megawatts,  the  grid  row  and  column  (location),  and  the 
availability  (decimal)  of  the  particular  generating  unit,  and  then  asks  for 
any  changes. 

In  the  second  section  of  generating  unit  data,  the  program  prints  the 
incremental  heat  rate  coefficients,  the  fuel  type,  the  fuel  cost  in  dollars 
per  million  BTU,  and  the  sulfur  content  in  percent.  (Note  the  incremental 

heat  rate  curve  for  each  unit  is  stored  as  a  quadratic  function  of  the  form 

2 
a+  P+yP  .)  After  printing  this  data,  the  program  asks  for  any  changes. 

In  the  third  and  final  section  of  generating  unit  data,  the  program 
prints  the  stack  height  in  meters,  the  stack  radius  at  the  top  in  meters, 
the  stack  gas  exit  velocity  in  meters  psr  second,  the  fuel  heating  value 
in  BTU's  per  pound,  and  th«  minimum  sulfur  content  in  percent.  The  program 
next  asks  for  any  changes. 

This  completes  the  generating  unit  data  for  the  particular  unit  the  user 
specified.  The  progra.Ti  aqain  asks  the  user  it  he  wishes  to  change  any  of 
the  generating  unit  data.  If  so,  the  process  described  above  is  repeated. 
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This  second  category  of  input  (generatiny  unit  data)  is  continued  until 
the  user's  response  indicates  he  no  longer  wishes  to  modify  any  generating 
unit  data. 

The  input  data  are  summarized  in  Table  4,  v/here  the  dotted  lines 
separate  the  sections  described  above.  The  ability  to  change  or  modify 
any  of  the  input  data  shown  in  Table  4  makes  this  program  a  \'ery   powerful 
tool  for  evaluating  operating  strategies  and  for  performing  planning 
studies  and  simulations. 

The  reader  is  referred  to  the  User's  Guide  to  ENVECO  for  a  detailed 
description  of  the  output.  We  will  only  mention  here  some  of  the  more  im- 
portant aspects  of  the  output  from  the  program.  We  will  end  this  section 
with  a  partial  list  of  some  of  the  information  the  program  can  provide: 

1)  The  results  of  the  effective  stack  height  computations 
(Chapter  3)  for  those  units  burning  coal  or  oil. 

2)  The  basis  for  which  each  group  of  units  was  committed,  i.e., 
on  economics  or  air  quality. 

3)  The  commitment  order  of  the  generating  units  within  each  group. 

4)  Whether  or  not  the  block  (see  Figure  15)  "compute  new  sulfur 
contents,"  was  invoked,  and  if  so,  the  new  schedule  of  sulfur 
contents. 

.  5)  The  generating  unit  fuel  costs,  emissions,  and  emission  rates. 

6)  The  total  group  and  system  fuel  costs  and  emissions. 

7)  The  predicted  grid  concentrations  of  sulfur  dioxide. 

8)  The  24  hour  l-st  and  2-nd  highest  predicted  concentrations. 

9)  The  3  hour  l-st  and  2-nd  highest  predicted  concentrations. 
10)  The  additional  fuel  costs  (if  any)  for  committing  some  of 

the  units  on  the  basis  of  air  quality. 
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Category  1 
System  Data 


Table  4.  Input  Data  Summary 


Category  2 
Generating  Unit  Data 


LDC  array 


Wind  rose  frequencies 

Height  of  stable 
Layer  aloft 

Time  period 
Grid  size 
Number  of  rows 
Number  of  columns 

10-m  wind  speed 
Ambient  temperature 
Number  of  groups 
Number  of  segments 

Number  of  units  in  group  1 

Number  of  units  in  group  2 

Number  of  units  in  group  3 

Number  of  units  in  group  4 


Maximum  output 
Grid  row 
Grid  column 
Unit  availability 

IHR  coefficients: 

Alpha 

Beta 

Gamma 
Fuel  type 
Fuel  cost 
Sulfur  Content 

Stack  height 
Stack  radius 
Stack  gas  velocity 
Stack  gas  temperature 
Fuel  heating  value 
Minimum  sulfur  content 
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5 . 3  Conclusion  and  Future  Work 

We  have  presented  a  new  approach  for  reducing  ground  level  concen- 
trations of  sulfur  dioxide  due  to  electric  power  plants.  Since  the  approach 
involves  time  frames  on  the  order  of  a  month  to  a  year,  it  is  not  an 
intermittent  control  scheme,  and  therefore,  also  reduces  emissions  of  sulfur 
dioxide.  It  is  an  approach  to  be  used  as  a  supplement  to  and  not  in  lieu 
of,  constant  emission  controls.  We  have  concluded  that  this  approach,  as 
an  interim  {up  to  1985)  control  technique,  should  be  accepted  among  other  i 
approaches  in  keeping  with  the  best  interests  of  our  nation.  In  this  respect  f 
we  will  conserve  valuable  low  sulfur  fuels,  make  them  available  to  sources  | 
more  likely  to  violate  the  primary  standards,  and  meet  the  primary  and  | 

secondary  standards  in  a  reasonable  and  acceptable  manner;  until  such  time       I 
as  widespread  use  of  constant  emission  controls  are  possible,  both  from 
an  engineering  and  socio-economic  standpoint. 

A  novel  solution  to  the  gradient  transport  dispersion  equation  has  been 
presented  and  shown  to  lead  to  the  familiar  and  widely  used  Gaussian 
dispersion  model.  Hanna's  (1972)  equation  for  predicting  average  concen- 
trations near  the  source  has  been  corrected  in  this  work;  and  a  more 
fundamental  relationship  connecting  the  gradient  transport  theory  to  the 
statistical  theory  of  dispersion  has  been  presented.  The  extension  of  the 
theory  of  production  costing  to  include  certain  nonlinear  functions  presented 
in  this  work  has  enriched  the  method. 

Finally,  in  addition  to  the  usefulness  of  this  approach  in  evaluating 
operating  strategies,  it  will  likely  find  many  uses  as  a  planning  technique 
by  system  planners.  In  view  of  this,  suggestions  for  future  work  include 

the  extension  of  the  theory  of  production  costing  to  include  any  nonlinear       • 

i 

function  and  useful  techniques  for  considering  partial  outages  of  generating 
units. 


APPENDIX  A 

Since,  as  is  mentioned  in  Chapter  4,  concentrations  are  approximately 
lognormally  distributed,  we  must  have  an  understanding  of  some  of  the 
properties  of  the  lognormal  distribution  that  are  used  in  Section  4.4. 
In  the  interest  of  completeness,  these  properties  are  discussed  in  this 
Appendix.  Also,  while  the  lognormal  distribution  is  not  uncommon,  the 
reader  is  probably  not  as  familiar  with  its  properties  as  he  is  with 
the  properties  of,  say,  the  normal  distribution. 

The  lognormal  distribution  has  had  applications  in  many  areas  of 
engineering  (see  Aitchison  and  Brown  (1957)).  One  use  of  the  lognormal 
distribution  has  been  for  problems  whose  variate  is  limited  to  the  positive 
real  axis.  A  function  that  is  normally  distributed  after  the  transformation 

y  =  In  X  (A-1) 

where  In  is  the  natural  logarithm,  is  said  to  be  lognormally  distributed. 
Its  probability  density  function,  f(y),  isgiven  as 

F(y)  ^  7^  exp  [-  y  i^fl 

t 
Its  cumulative  distribution  function  will  be 


F(y)  = 


ry 


f(y)dy 


Using  the  transformation  (A-1),  the  cumulative  distribution  becomes 


t 


In  Chapter  4  the  cumulative  distribution  function,  F'(-)>  had  a  nrime  over 
it  to  distinguish  it  from  the  backward  cumulative  distribution  function, 
F(-).  In  this  appendix  the  prime  is  unnecessary  and  so  we  drop  it. 
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F(ln  x)  =  f  "  "^  -^  exp  [-  \   (ln-^-^Ji.)2jdx  (A-2) 

JO    X/2TO       "     "^ 

Equation  (A~2)  can  be  rev/ritten  as 

/•In  X 
F(ln  x)  =  I    f(x)dx 
^0 

where 

F(x)  =  -1^  exp  [-  i(IlLx_:iJi)2]  (A-s) 

X/2TTa 

Equation  (A-3),  then,  is  the  probability  density  function  for  a  lognormal 
distribution.  It  is  plotted  in  Figure  15  for  i.i=0  and  a=l . 

Equation  (A-3)  is  seen  to  be  characterized  by  two  parameters, 
y  and  a.  We  need  estimators  of  these  parameters  to  be  able  to  describe 
a  lognormal  distribution.  We  will  use  the  maximum  likelihood  estimator. 
For  the  estimator  of  y  we  have 

L(y)  =  f(x-|,X2,...,x^;y)  (A-4) 

where  L(*)  is  the  likelihood  function,  x,  ,X2,. . .  ,Xj^  are  sample  values 
from  a  sample  size  n,  and  f  (x-,  ,X2» . . .  ,x  ;y)  is  the  joint  density  function 
of  the  n  samples  and  the  parameter  to  be  estimated,  y.  The  maximum 
likelihood  estimator  assumes  statistical  independence  and  so  equation 
(A-4)  becomes 

L(y)  =  f(x^;y)f(x2;y)...f(x^;y) 
where 

1        r   1  ."'^  'i-^-,2. 


f(xr'Vi)  =  — iz-  exp  [■-  ^  ( — ) 

1 


..,^;         ..p  u-  2  _   ^ 


'^^"  ■        1     n   T       -^  In  x.-y  , 
L(y)  =   (~~r)   n  (-L  exp  [-  \  (— -M^]) 


V 


li'O         1=1   1 


»mJa.c,^^S£. 
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To  determine  the  maximum  likelihood,  it  is  necessary  to  take  the 
derivative  of  L  with  respect  to  y  and  set  the  result  to  zero.  However,  the 
same  value  of  u  that  maximizes  L  will  maximize  In  L.  The  use  of  the 
natural  logarithm  greatly  simplifies  the  computation.  Thus 

n        1  n  In  x.-y  2 

In  L(y)  =  -n  In  (/2^a)  -  I     In  x  -  ^  1   ( ^) 

i=l    ^   ^  i=l    ° 

Then,  taking  the  derivative  and  setting  it  equal  to  zero,  we  have 

.  -,  n  In  x.-y 

£  [in  L(y)]  =  1 ;  (-^)  =  0 


from  which 

n 

I     In  x.  =  ny 


i=l 
or 


1  Y  In  x.  (A-5) 


u  =  -  y  In  X. 
n  i^l 


Notice  that  the  maximum  likelihood  estimator  of  y  is  seen  to  be  analogous 
to  the  sample  mean,  except  that  the  logarithm  of  the  sample  is  used. 

If  exponentials  of  both  sides  of  equation  (A-5)  are  taken,  we  have 

1  n 
exp(y)  =  exp  {j-    I     In  x.) 
"  i=l    ^ 

1    " 
exp(y)  =  exp  [-  ln(  n  x.)] 


i=l 
n 


exp(y)  =  [  n  x.]^/"  (A-5) 

i=l 

The  right  hand  side  of  equation  (A-6)  is,  by  definition,  the  geometric 

mean,  m  ,  thus 

m  =  exp(y) 


or 


y  =  In  m  (A-7) 

9 

Thus,  equation  (A-3),  the  lognormal  probability  density  function,  becomes 


HI 


T  -,     In  x.-ln  m     ^ 


x/2TTa 

Similarly  for  the  estimator  of  o,  the  likelihood  function  is 

,  n      -,  1     In  x.-ln  m     ^ 

L(a)   =   {-^f     n    f  exp[-l( V-^)2] 

^^a       i-1   X.  ^  u 

Taking  the  logarithm  yields 

n        ,  n  In  x.-ln  m  ^ 

In  L(a)  =  -n  ln(/2^a)  -  J  (In  x  )-  i  I  ( '- ^f 

i=l    ^   "^  i=l     "^ 

and  differentiating  and  setting  the  result  to  zero,  we  have 

,                n  In  x.-ln  m  ^ 
^[1nL(a)]  --^-     I   ( '-^ ^)^  =  0 

or 

n  ,ln  x.-ln  m  .2     1/2 

a  =  [  I  ^ ^^^ 2_)  ]  (A-9) 

i  =  l 

Similar  to  the  relationship  between  the  geometric  mean  and  y 

m  =  exp(u) 
9 

we  define  the  geometric  standard  deviation  as 

s„  -   exp(a)  (A-10) 

9 

where  we  have  chosen  s  rather  than  a  to  be  consistent  with  Larsen  (1959) 

g        g 

Equations   (A-9)  and   (A-10)  together  relate  the  standard  geometric 
deviation  to  the  samples.  Note  the  similarity  of  equation   (A-9)  to 
the  sample  standard  deviation  of  a  normal  function.  Similar  to  equation 
(Ar-7)   we  rewrite  (A-10)  as 

a  =  In  s  (A-11) 

g 
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and 


In  m. 


(A-7) 


where  equation  (A-7)  is  repeated  here  for  convenience. 

The  need  for  the  geometric  mean  and  standard  deviation  becomes 
clearer  in  Section  4.4.  To  complete  this  appendix  ,  v/e  need  to  determine 
the  mean,  median,  and  mode  of  the  lognormal  distribution.  Recall  that 
for  the  normal  distribution,  they  are  all  three  the  same  and  equal  to  y. 

The  mean,  m,  or  expected  value,  E[x],  is 


m=E[x]  = 


E[x] 


0 


xf (x)dx 


1  1   „   r  1  /■lnx-iJ\2-,  , 

X exp  [-  p-  ( -^  ]  dx 


Consider  the  change  of  variable 


y  =  In  x 


Then 


1 


E[x]  =  E[exp{y)]  =    exp(y)  -33-  exp  l-  2^  ^ 

J  -co      /Z-rra 


[-  l(^)^]dy 


which  can  be  rewritten  after  grouping  of  some  terms  as 


E[x]  =  -—  ["  exp  [-  X.   (y2-2(y+a''^)y4-y2)]dy 


2a 


Now,  if  we  complete  the  square  of  the  terms  in  parentheses,  and  let 


a  =  \x-^Q 


we  get 


•03 
E[x]  =  expdi-f  ^f^)     -^-  I       exp[-  k^^-)^-]dy 


2Tro  ^  -°o 


i.      c 


!43 


or 


1  2 
E[x]  =  exp  (y+  2-  a  ) 


(A-12) 


since  the  term  in  brackets  evaluates  to  one. 

The  median,  by  definition,  is  the  50th  percentile  position.  That 
is,  the  abscissa  value  for  which  the  cumulative  distribution  function 
equals  1/2.  From  equation   (A-12)  this  becomes 
In  X 


1 


'      ^      exp  [-  i  (ll^J-=Ii)2]dx 


2tjo 


So  we  want  to  find  the  abscissa  value,  In  x  ,  for  which  the  right  hand 
side  evaluates  to  1/2.  Again,  we  use  the  logarithm  transformation 


so 


and 


y  =  In'  x 


y.  =  In  X 
•'o      0 


/o 


y2?( 


1      r  l/y-y>2-, 
—  exp  [-  ji^""-)   ] 


■no 


2'  o 


dy 


which  is  in  the  form  of  a  normal  distribution.  For  a  normal  distribution, 
we  know  the  median  is  y.  Thus 


Yq  =  M 


from  which 


Xq  =  expfp) 


lie  now  corns  to  one  of  the  most  important  properties  of  a  Icgnormally 
distributed  function,  because  from  equation  (A.-7) 


p  =  in  m 
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or 

g       r'-^/ 

and  therefore  the  median  equals  the  geometric  mean,  lie  will  express 
this  as  an  equation  since  we  will  need  to  refer  back  to  it.  Thus 

Median  =  x  -  m  (A.-13) 

0   9 

Finally,  the  mode  is  the  maximum  value  of  the  probability  density 
function  (our  function  is  unimodal ) .  Differentiating  equation  (A-3) 
and  setting  it  equal  to  zero  yields 

VZ-na 

or 

_T  -  (In  x)-p  ^  Q 

a 

2 

(In  x)-y  =  -a 

In  X  =  u-a 

2 
X  =  exp(y-a  )  (A-14) 

Equation  (A-14)  is  the  equation  for  the  mode  of  a  lognormal  distribution, 
Figure  16  shoves  a  lognormal  density  function  for  y^O  and  0=1,   for  which 

mode  =  exp  (-1 )  =  .368 
median  =  exp(O)  =1.00 
and 

mean  =  exp  (1/2)  =^  1 .65 
and  these  values  are  shov;n  in  the  figure. 


145 


X3 

-Q 
O 

s- 

D- 


'^o'^^  median  "^^a" 


Concentrations  In  Micrograms  per  Cubic  Meter 


Figure  16.  Probability  Density  Function  for  a  Lognormal 
Distribution. 
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